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INTRODUCTION 



I . 

The trend in packaging of electronic components is to 
place more power in an ever decreasing space. Computers that 
occupied an entire room 20 years ago, now occupy a small 
corner of that room. With the size of chips decreasing, the 
volumetric power generation (q'") is of the order of MW/m 3 in 
magnitude . 

This can lead to high temperature operation of the chips 
with a resulting degradation in performance and or failure. 
The military has set the maximum temperature of a 
semiconductor junction to between 100 and 110 ’C [Ref. 1], 

This investigation deals with the steady state temperature 
distribution and resulting thermally induced stresses in a 
tri-material electronic package as shown in Figure 1.1, 
dimensions not shown to scale. The package consists of a 

semi-conductor (chip) attached to a substrate material by a 
solder joint. The chip is a source of thermal energy. The 
package is exposed to a convective environment on all sides. 



\ 


1 

CHIP 


1 "1 





SOLDER 


: 


SUBSTRATE 

1 

1 



Figure 1.1 Tri-material electronic package 
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With the increase in volumetric power generation the 
overall temperature of the package is increased. The stresses 
in the package arise from two sources, (a) thermal strains 
( €j ) caused by temperature gradients within a material and (b) 
differences in the materials mechanical properties, 
coefficient of thermal expansion (a) and stiffness (E). At 
material interfaces, displacement continuity requires 
expansion or contraction accommodations between materials. 
There are two such interfaces where this will occur, at the 
chip/solder interface, and at the solder/substrate interface. 

A . BACKGROUND 

Thermal stresses of bi-material assemblies were 
investigated by Timoshenko in a 1925 paper [Ref. 2], where a 
bi-metallic thermostat assembly was subjected to a uniform 
temperature field. Timoshenko obtained analytical expressions 
for the bending stresses through the cross section of the 
bimetallic strip. He also noted that the 'distribution of 
shearing stresses along the bearing surface cannot be 
determined in an elementary way, that they are of a local type 
concentrated near the ends of the strip', and that the 
shearing stresses can be of the same magnitude of the bending 
stresses. He mentions, in passing, the existence of 'local' 
normal stresses between the assembly interfaces but says 
nothing of their magnitude. Goland and Reissner [Ref. 3] 
determined the shear and peeling stresses at the interface 
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of a cemented lap joint, where the peeling stress tends to 
pull the materials apart. 

In a comprehensive review of the subject of tri-material 
assemblies, Suhir [Ref. 4] took the analysis one step further 
by investigating the stresses in a tri-material system due to 
a thermal environment. His model is based on an assembly 
fabricated at elevated temperatures and subsequently cooled. 
Specific restrictions were placed on the tri-material assembly 
in a uniform temperature field. He obtained analytical 
solutions of the problem. For the present study the general 
tri-material problem is solved numerically. In particular the 
tri-material system is that of an electronic package. 

B. PROBLEM DEFINITION 

The problem was partitioned into two parts, a thermal 
analysis, whereby the heat conduction equation was solved, and 
a thermal stress analysis whereby the stresses produced by the 
thermal field were obtained. Both problems were solved 
numerically by the Finite Element Method (FEM). For the 
numerical analysis, symmetry is invoked along the centerline 
of the system as shown in Figure 1.1. The steady state 
thermal study focuses on the effects material properties, 
geometric configuration, convective cooling, and encapsulation 
have on temperature distribution in an electronic package. 
Once the temperature field is known the data is input into a 
program that assigns temperatures to another FEM mesh which is 
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used in the stress analysis. The stresses examined include 
bending stress at the centerline of the package (a b ), the 
bending (a bu ) , normal (o 0u ), and shear (r u ), stresses along the 
chip/solder interface, and the bending (a bL ), normal (a oL ), and 
shear (r t ) stresses at the solder/substrate interface. 
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II. FEM DEVELOPMENT OF THERMAL PROBLEM 



A. PROBLEM STATEMENT AND ASSUMPTIONS 

As noted in the introduction, the purposes of this 
investigation are to (a) determine the variation of 
temperature in the electronic package and (b) determine the 
stresses associated with these temperatures. The results here 
are the steady state results. 

The steady state thermal study focuses on effects that 
certain key parameters have on system behavior. They include, 
the effects total power (P) in the chip, material 

combinations, convective cooling, geometric configuration and 
encapsulation have on the temperature distribution in the 
electronic system. 

Finite element codes were used to solve for the 
temperature distribution and stresses of the electronic 
package. In order to limit the mathematical complexity of the 
problem the following approximations and specific assumptions 
for FEM program HEATSTEADY in Appendix A. are made as follows: 

• The thermal contact resistances due to 'imperfect' contact 
at the solder interfaces are negligible. 

• The thermal conductivity and coefficient of thermal 
expansion (k and a) are considered constant over the range 
of temperatures encountered in the study. 
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• Because the temperatures are not very large radiation heat 
transfer between the package and the surrounding enclosure 
is negligible. 

• The convection heat transfer coefficient (h) is constant 
over the entire edge upon which it acts . 

All of these assumptions are reasonable within the scope 
of this study. Once the temperature distribution of the 
package is obtained the position and temperature of all the 
Global nodal points are entered into program Filter in 
Appendix B. to obtain temperatures that correspond to FEM 
program Weld in Appendix C. to determine the stresses 
developed in the package. The development of program WELD is 
discussed in Chapter III. 

B. FEM FORMULATION OF THERMAL PROBLEM 

The Galerkin FEM is an approximation method which 
transforms a linear partial differential equation into a 
system of linear algebraic equations. Using the two- 
dimensional heat equation [Ref. 5]: 



V • (£V T) + q'" = 0 



( 2 . 1 ) 



and boundary conditions 



Cauchy : -k^=h(T-Tj 

on 



( 2 . 2 ) 



Nuemann : 




(2.3) 
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and thus 



W-(kVT) = kV^T = k 



&T + yr 



( 2 . 4 ) 



dx 2 dy 



2 



where k is the thermal conductivity, T is the temperature, and 



distribution can be determined at system nodal points of the 
electronic package. For this method, a linear triangular 
shape function (N f ) which possesses the Kronecker Delta 
property (equation 2.5): 



is used. The linear shape functions maintain function 
continuity throughout the domain. 

Heat flux continuity at material interfaces is built into 
the FEM formulation. The thermal FEM grid and boundary 
conditions are given in Figure 2 . 1 and Figure 2 . 2 
respectively . 



q'" is the volumetric heat generation. The temperature 



N £ (. nodej ) = 6^ if i = j then = 1 

if i * j then N± = 0 



( 2 . 5 ) 
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FEM thermal grid 
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Figure 2.2 FEM thermal problem boundary conditions 



9 



An approximate solution, t, for temperature, T(x,y), is 



formed as follows: 



T ~ t = {M T (f} 



( 2 . 6 ) 



where T is the exact solution of the heat equation in 
continuous space, t is the approximate solution in discrete 
space, {N } T is the transpose of a column vector of the 
triangular linear shape functions, and {T} is the vector of 
nodal temperatures . 

After the approximation is formulated, the next step is to 
form the Residual, R, as follows: 



where q'" is the heat generation term and $6 denotes the 
differential operator which in the case of the heat equation 
is defined by: 



R = S£(t) - g'"{s) 



( 2 . 7 ) 



S£ ( t) = V-(kVt) 



( 2 . 8 ) 



With this substitution, the residual becomes: 



R = [V-(kV - q'"{s) 



( 2 . 9 ) 



From the residual, the Galerkin Equations are formed: 




( 2 . 10 ) 



10 



where D denotes surface integration and where {0} is the null 
vector. Further substitution for R into the Galerkin vector 
equation results in: 

f In) [V-(kV ({n) t {t)) ) ds - f { N)q'"(s)ds = (o) (2.11) 

J D J D 

To solve the Galerkin Equation, Green's Theorem is evoked 
which yields: 

f In} (kVT) dB - f [V{at> (icV ( {n) t {t}) ) ] ds + [ {N)q'"ds= (0> 

where B subscript on the integral denotes evaluation of these 
integrals around the boundary of domain D of the electronic 
package. When integrated the boundary integral and the heat 
generation integral form the "force" vector {F}. The middle 
integral forms the [A] matrix. Thus the heat conduction 
differential equation becomes, 

[A]{T) = {F] (2.13) 

where {T} contains the vector of nodal temperatures. 

C. THERMAL PROGRAM DESCRIPTION AND VALIDATION 

The main thermal program, HEAT. FOR, a VAX Fortran 77 Code 
for the FEM thermal analysis is contained in Appendix A. The 
program begins by reading in a data file DATA3.DAT. This 
input data is then processed through several subroutines that 
perform the FEM analysis. 
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1. Subroutine Input Description 

The following information is supplied to the Input 



subroutine 


from a data file. 


Input file 


form: NC0L1 ,NC0L2 ,NROW, 

XPOS , YPOS, DIVX, DIVY, 

THERCON , GEN , QTPR , 

DIVSUB, DIVXAIR, DIVYAIR, 

NBCLOW , NBCUP , NBCVERT , NBCSUB , TAMB , 
HLOW , HUP , HVERT , H SUB 


Parameter 


Description 


NC0L1 


The number of columns of substrate in the FEM grid. 


NC0L2 


The number of columns of chip in the FEM grid. 


NROW 


The number of rows in the FEM grid. 


XPOS 


Starting position for new grid density along X 
axis . 


YPOS 


Starting position for new grid density along Y 
axis . 


DIVX 


Spacing of grid density along X axis. 


DIVY 


Spacing of grid density along Y axis. 


THERCON 


Thermal conductivity (k). 


GEN 


0 = no heat generation term 

1 = heat generation term 


QTPR 


Amount of volumetric heat generation (cf") . 


DIVSUB 


Number of rows in FEM grid composed of substrate. 


DIVXAIR 


Number of columns in FEM grid composed of air. 


DIVYAIR 


Number of rows in FEM grid composed of air. 


NBCLOW 


0 = insulated boundary at the bottom of substrate. 

1 = convective boundary at the bottom of substrate. 


NBCUP 


0 = insulated boundary at the top of chip. 

1 = convective boundary at the top of chip. 
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NBCVERT 0 = insulated boundary on the side of the chip. 

1 = convective boundary on the side of the chip. 



NBCSUB 0 = insulated boundary on top of the substrate. 

1 = convective boundary on top of the substrate. 



HSUB 



HVERT 



HUP 



HLOW 



Convection coefficient corresponding to NBCLOW. 
Convection coefficient corresponding to NBCUP. 
Convection coefficient corresponding to NBCVERT. 
Convection coefficient corresponding to NBCSUB. 



2 . Subroutine Grid Description 

Subroutine GRID takes the input information and 
constructs a grid of horizontal and vertical lines on a 
geometry as shown in Figure 2.1. This grid will be used to 
form the triangular elements used in the thermal FEM analysis. 
This subroutine allows the user to generate a variety of 
meshes and/or refinements and permits validation of grid 
independence very quickly. 

To optimize computer time the subroutine is able to 
generate different grid meshes with a variety of spacing 
options as shown in Figure 2.3. Fine meshes are used at 
convective boundaries and material interfaces. A coarse mesh 
is used within a material away from boundaries. The 
subroutine allows the user to independently vary the mesh 
density both horizontally and vertically for up to nine 
different densities in either direction. The number of mesh 
densities can be increased by simply adding a number of "IF 
THEN" statements at the top of the program (line 34). 
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fine 



coarse 



coarse 





fine 

Figure 2.3 Fine and coarse vertical mesh and fine and coarse 
horizontal mesh 

Local nodal points are assigned to each element, one 
at each corner, counterclockwise as shown in Figure 2.4. 




Figure 2.4 Triangular element numbering 

Next the subroutine generates a correspondence table 
which identifies the correspondence between local nodal points 
and global nodal points. Next, each triangular element is 
identified with a specific material and assigned a thermal 
conductivity (k) and whether or not it will have a volumetric 
heat generation term (q'"). Lastly boundaries are assigned 
according to Figure 2.2. 
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3. Subroutine Heatmat, Heatmod and Lsarg Description 

Subroutine HEATMAT constructs the global matrix [A] 
and the heat generation vector {F} according to the 
correspondence table made in subroutine GRID. 

Subroutine HEATMOD is not used in this program but can be 
used to modify the global matrix [A] for a Dirichlet boundary. 
HEATMOD replaces the Galerkin equation for the specific nodal 
points with the specified temperatures. LSARG is an equation 
solver provided by the IMSL Library used to solve equation 
2 . 13 . 

4. Subroutine Output Description 

Subroutine OUTPUT creates data files and performs an 
energy balance between the thermal energy which is convected 
at the boundaries and the thermal energy generated by the chip 
according to: 

E hi < r t - T.) L t = q"'A cU , (2-10 
2=1 



where h A and 1 A correspond to the convection coefficient and 
length of each convective boundary, T t corresponds to the 
average temperature across the boundary and A chip is the area 
of the chip. This heat balance is conducted for each case 
study. 
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D. VALIDATION I. DIRICHLET AND CADCHY BOUNDARY 



A two dimensional body subjected to both Dirichlet and 
Cauchy boundary conditions was taken from [Ref. 6] and run on 
program HEAT to validate the programs ability to handle both 
boundary conditions. The results of program HEAT were 
identical to Reference 6. 

E. VALIDATION II. GRID INDEPENDENCE 

In order to ensure that the numerical errors are kept to 
a minimum the size of the mesh was decreased with incremental 
steps until there was a less than 1% change in the solution 
from one grid to the next finest grid (grid independence) . 

The final mesh chosen has 125 degrees of freedom (DOF) 
with 202 elements. The mesh which provided less than 1% 
difference had 161 DOF with 268 elements. The total power of 
the chip used is P = 40 W. The total power convected for the 
DOF = 125 grid was P conv = 39.95 W. The total power convected 
for the DOF = 161 grid was P conv = 39.94 W. The percent 
difference in total power P from the coarse mesh to the fine 
mesh is 0.03%. 

F. VALIDATION III. ENERGY BALANCE 

As mentioned previously, there is an energy balance in 
subroutine OUTPUT which is given by equation 2.14. This 
checks program HEAT's ability to correctly handle a generation 
term. The total power for the majority of case studies was 



16 



equal to P = 1 W. The average convective power output was 



q conv = .99785 W for a percent difference of 0.21%. 
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III. FEM STRESS DEVELOPMENT 



The finite element stress formulation utilizes a recently 
developed element which provides for axial and lateral 
displacement continuity. The results of the stress code are in 
good agreement with existing solutions [Ref. 2,4] to problems 
with uniform temperature fields. 

A. FORMULATION OF THE MODEL 

A brief description of the FEM formulation for stresses 
follow. In Figure 3.1, each element has six degrees of 
freedom; axial displacements at the four corner points, and 
lateral displacements at the two ends. An advantage of the 
element is that axial and lateral displacement continuity 
results . 




Figure 3.1 A typical element with 6 degrees of freedom 
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The axial displacement field u(x,y) is assumed to be linear in 
both the axial and transverse directions. That is, 

2 

u (x, y) = '£N i (x)[H 1 (y) u? +H 2 (y) uf] (3.1) 

i=l 



Superscripts b and t on the nodal displacements refer to 
bottom and top displacements respectively. The linear shape 
functions N i and H i are given by 

H x (y) = h 2 = Z n 2 (x) = n 2 = z 



The lateral displacement field is given by 

2 

v(x) =52 N i (x) V i 

i = 1 



(3.2) 



The strain-displacement relations are, 



du (y) r b ts\ H 2 (y) r t ti ,, 

e * = —h - u iJ + —j— [ u * - ui] < 3 - 3 > 



and 



N 1 (x)r t M N 2 (x) r c „ Vj ~ V 

Vxy = K - U l] + [ U 2 - U 2\ + — f - < 3 4 > 
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Defining the axial displacement vector as 



= { uf Ilf u 2 u 2 c ) (3.5) 

and the stiffness matrix as 

[iC B ] = f^f^EiBFdydx (3.6) 



where the B vector is, 



, , dN ' dN. 

{B) T = ( —iH. 



dx 



dx 



dN, 

dx 



H, 



dN, x 



(3.7) 



and the force vector due to temperature as, 

(NJ = ^ f 1 [ h iB}EaATdydx (3.8) 

611 o Jo 

gives the bending matrix eguations as, 

[jg {6^} = (F^} (3.9) 

Eguation (3.9) defines the bending behavior. Behavior due to 
shear is obtained as follows. Define the row vector of 
displacement degrees of freedom as 

{6j T «(u? uf v ! u, u, v 2 ) (3.10) 

The shear stiffness matrix is given by 
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[K s ] = f^f^B^GiBOdydx 



(3.11) 



where 



(BO = (N 2 



dH 1 

dy 




dN, dH, 



dbl* l . _ — _ . 

N 2 ^ -^) (3.12) 

z dy dx 



which gives the equations for shear behavior 

[K 3 ] (6 s } = (0} (3.13) 

The matrix equations for bending and shear behavior are 
combined to give the stiffness equations for the system. 

B . PRELIMINARY RESULTS 

Assessment of the present model was obtained from 
comparisons with previous analytical models [Ref. 2,4]. Exact 
agreement for bending stresses was obtained when the present 
model was applied to a sample bimetallic case solved in 
[Ref. 2]. For a tri-metallic case the present model produced 
shear stress results in reasonably good agreement with results 
obtained from [Ref. 4]. 
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IV. STEADY STATE TEMPERATURE ANALYSIS 



A. OVERVIEW OF STEADY STATE TEMPERATURE ANALYSIS 

The analysis begins with an FEM approximation of the 
temperature distribution of the tri-material electronic 
package. First the chip volumetric power generation is 
examined to determine the effect an increase in q'" has on the 
chip temperature. Keeping the power of the chip and the level 
of convective cooling constant, several packaging materials 
are combined to determine how the temperature distribution 
varies with the variation in the material thermal 
conductivities. This is followed by a convection cooling 
analysis, in which temperatures and temperature gradients are 
examined for a common electronic package while the convection 
coefficient h varies. The package geometry is then studied, 
holding the overall package height equal and varying the chip 
size as a percentage of overall package height. This is done 
for the case of constant heat flux, and for the case of 
constant volumetric heat generation. Lastly the effects of 
placing a protective coating "encapsulation" on the chip is 
studied to determine the variation in temperature distribution 
from a chip which has not been encapsulated. 
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B. EFFECT OF VOLUMETRIC POWER GENERATION ON TEMPERATURE OF 



THE CHIP 

The heat fluxes of electronic packages is expected to 
reach nearly 10 6 W/m 2 by the year 2000. This is two magnitudes 
in power over the chip studied here. The total power output 
(P) of the chip studied is one watt. This translates to a 
volumetric power generation (q'") of 20 MW/m 3 for a device 5mm 
x 5mm x 2mm in dimension or a heat flux (q") of 4xl0 4 W/m 2 . 
The dimensions chosen are of a typical electronic device size. 

The temperature results throughout the study can be 
modified for any power output due to the linearity of the 
generation term in the heat equation. In Figure 4.1 the 
temperature of the chip is plotted against chip volumetric 
power output (q'''). All other quantities must be constant 
for this to remain true, including boundary conditions and 
other materials within the package. For the particular 
geometry studied, temperature (T) as a function of volumetric 
power output q'" W/m 3 is given by equation 4.1 

T = m * q'" + T m ( 4 . 1 ) 

where the change in temperature per unit change in volumetric 
power output is b = 4.65 'C/MW/m 3 and T. is the ambient 
temperature of the convective fluid. 
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120 




Figure 4.1 Chip temperature as a function of volumetric power 
C. EFFECT OF VARIOUS MATERIAL COMBINATIONS 

The materials used in this investigation along with their 
key thermal properties are listed in Table 4.1. The set of 
materials that comprise the electronic tri-material package 
focused on for the majority of this study are the silicon 
chip, a Pb-Sn solder, along with several substrate materials 
including, Epoxy Fiberglass, Polyimide Fiberglass, Alumina, 
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and Aluminum Nitride. These are the most common materials 
used in electronic packaging. The other chip and solder 



materials 


considered in this 


study account for large 


variations 


in thermal conductivity and/or 


thermal expansion 


coefficient (a). 






TABLE 4.1 


ELECTRONIC PACKAGE HATERIALS AND 


PROPERTIES 


Use 


material 


k (VI /mK) 


a (10' 7 /°C) 


Chip 


Silicon 


150 


27 




Gallium Arsenide 


58 


57 


Substrate 


Epoxy Fiberglass 


0.16 


140 - 180 




Polyimide Fiberglass 


0.35 


120 - 160 




Alumina 


18 


60 




Aluminum Nitride 


230 


33 


Solder 


Pb(5%)/Sn 


63 


290 




Au( 20% ) /Sn 


57 


159 


Air Gap 




0.026 


— 



Throughout the study within this section the following 
parameters are held constant; volumetric power (q'") of the 
chip, ambient temperature T„ of the convective fluid, and 
coefficient of convection cooling (h). Figure 4.2 shows the 
non-dimensional temperature profiles of the package at a non- 
dimensional distance, <£ x = 0.564. The plotted values of T n 
were obtained at the FEM grid nodes. The temperatures which 
lie along <t> K = .564 were chosen because they pass through the 
solder joint of the package. Thermal stresses will result not 
only from temperature gradients within a material, but also 
from the direct contact of the three different materials 
within the package which have large differences in their 
coefficients of thermal expansion (a). In Figure 4.2 the 
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abscissa is a normalized temperature T n given by equation 4.2, 



T n = 




(4.2) 



where T. is the FEM grid temperature. The ordinate is the 
non-dimensional y location cf> y that is given by equation 4.3. 

4> y = jj (4.3) 

where H is the total height of the package and y is the FEM 
grid global nodal location. 




NON DIMENSIONAL TEMPERATURE T * 111 * 

n T_ 



Figure 4.2 Non-dimensional temperature for silicon chip, Pb-Sn 
solder and various substrates. 
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All sixteen combinations of materials were analyzed for a 
temperature profile. There was no significant difference in 
temperature profile with regard to which chip or solder 
material was used in the package. For the volumetric power 
level used in this analysis the difference in thermal 
conductivity of the Silicon chip (k=150 W/m'C) or Ga/As chip 
(k=58 W/m°C) does not influence the temperature profile. A 
study was conducted in which the volumetric power was raised 
by four times the initial amount. This failed to produce any 
significant difference in temperature profile or magnitude. 

The temperature profiles show that the substrate is the 
governing material which determines the temperature profile 
a given electronic package takes on. The following symbology 
will be used in further analysis. A subscript of 1 refers to 
the chip, subscript 2 refers to the solder, and subscript 3 
refers to the substrate. Low thermal conductivity substrates 
used in this analysis have ratios of k 3 /k., < 0.0025. High 
thermal conductivity substrates used in this analysis have a 
ratio range of 0.10 < kj/k., < 1.6. 

A characteristic of a high k 3 /k., package is a uniform 
temperature field throughout all materials in the package. 
All materials within the package can conduct heat very well. 
With the low k 3 /k., substrates a temperature gradient is 
developed within the substrate in both the X and Y directions. 
As a result of not being able to conduct the heat effectively 
through the low conductive substrate, the chip temperature is 
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increased by approximately 40% over a package with a highly 
conductive substrate. 

The solders have thermal conductivities which are close to 
the chip thermal conductivity (k 2 /k, = 0.446) and do not effect 
the temperature profile. However, the coefficient of thermal 
expansion a between the two solders chosen vary by 55% thereby 
influencing the thermal strain profile significantly. The 
solder joint takes on the chip temperature in all cases 
because of its ability to conduct the heat generated by the 
chip. 

Thermal strain profiles are presented in Figures 4.3 and 
4.4 for the two different solders. In these figures, the 




THERMAL STRAIN € * (T-T ) <t(x10 S m/m) 
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Figure 4.3 Thermal strain profile for silicon chip, Pb-Sn 
solder and various substrates. 
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NON DIMENSIONAL Y LOCATION <f> 



abscissa is the difference of the nodal temperature and 
initial temperature multiplied by the coefficient of thermal 
expansion a given by 

e T = (T i - TJ * a (4.4) 
and the ordinate is given by equation 4.3. 




Figure 4.4 Thermal strain profile for silicon chip, Au-Sn 
solder and various substrates. 
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1. Low Thermal Conductivity Substrate Packages 



Further analysis is conducted by dividing the packages 
into two groups, low thermal conductivity substrates and high 
thermal conductivity substrates. The two substrates chosen, 
Epoxy Fiberglass (k 3 /k n = 0.00106) and Polyimide Fiberglass 
(kj/k., = 0.0023), are part of a group of substrates that have 
k 3 ranging from 0.12 W/nTC to 0.35 W/nTC [Ref. 1]. The 
temperature achieved throughout the chip for each of the case 
studies, Figure 4.2, are approximately constant, that is, not 
dependent on position. The variation in temperature of the 
chip was limited to less than 1.2% over the entire cross 
section. The solder achieved the same temperature as the chip 
due to the solders high thermal conductivity. The variation 
between the chip temperature and the solder temperature is 
less then 0.1%. 

Within the substrate temperature gradients 9, and © 

x y 

are developed in both the X and Y directions respectively. 
Temperature gradients are given by: 



dT _ q _ ( 1 ^ 

dx x (x i - x i+1 ) 



( 4 . 5 ) 



dT _ q _ ( Tj \ 

dy y (yj - y>i) 



( 4 . 6 ) 



where T^ and T- are FEM Global nodal temperatures, x i and y- 
are FEM Global nodal distances in the X and Y directions 
respectively . 
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These temperature gradients are only found in low 
thermal conductivity substrates and an explanation for this 
can be found by examining Fourier's Law, equation 4.7, 

q " = ~ k % = *- fc6 ' <4 - 7> 

where q" is the heat flux in the X direction. 

The solder acts as a sort of 'heat funnel' for the 
chip to transfer its thermal energy to the substrate since the 
air pocket under the chip acts as an insulator. By examining 
Fourier's law there are two ways to increase the heat flux q". 

First by increasing the thermal conductivity (k) of 
the material, and secondly by increasing the temperature 
difference (or temperature potential). The convection heat 
transfer from the chip surface to the surrounding medium is 
unable to transfer all the thermal energy in the chip 
generated by the power source q"' and therefore results in the 
conduction of heat flux through the solder. When the thermal 
conductivity of the substrate is low (as in the Fiberglass 
substrates), the temperature difference through the substrate 
must become proportionally larger to account for the heat flux 
out of the chip. For this reason, when the substrate has a 
high thermal conductivity, the temperature potential need only 
be very small to conduct the same amount of thermal energy. 

In the X direction the gradient is largest at the 
upper right hand side of the solder/substrate interface. As 
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Figure 4.5 shows, the gradient 0 X diminishes by roughly 72% 
when only 25% into the substrate. The temperature gradient 




Figure 4.5 © x in epoxy fiberglass substrate 

has decreased by 88% when half way through the substrate. 
This behavior can be accounted for by the fact the thermal 
resistance of the substrate in conduction is greater then the 
thermal resistance of the convection surface at the top of the 
substrate. Equivalent thermal resistances in conduction and 
convection are given by: 
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< 4 . 8 ) 



R 



cond 



_AX 
k A 



where Ax is the distance between two points of interest 



R 



conv 



1 

h A 



( 4 . 9 ) 



and A is the area in which the heat flux is acting. 

Taking the area under the solder joint for conduction 
and the top of the substrate to the right of the solder joint 
for convection the following equivalent thermal resistances 
are obtained: 



R. 



cond 



0 . 002 

0.16 (0.00088) (1) 



R 



conv 



200 ( 0 . 001 ) ( 1 ) 



14.2 ° C/W 
5 °C/W 



( 4 . 10 ) 

( 4 . 11 ) 



which shows clearly that the thermal energy has approximately 
1/3 less resistance at the convective boundary versus 
conduction through the substrate. 

The gradient in the X direction on the left hand side 
of the solder joint is only 41% of that on the right, as shown 
in Figure 4.6. This again is a result of a higher thermal 
resistance in conduction within the substrate. With the 
insulated boundary condition imposed in the FEM formulation, 
the gradient in the X direction at X=0 and X=L is necessarily 
zero. 
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Y LOCATION THROUGH SUBSTRATE (mm) 



TOP OF SUBSTRATE 




Figure 4.6 0 X in epoxy fiberglass substrate 

The temperature gradient in the Y direction, 0 y , is 
largest at 49.02 'C/mm at the solder/substrate interface as 
shown in Figure 4.7. Away from this interface the gradient 
decreases to approximately 20 'C/mm and remains constant 
throughout the entire substrate. 
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Y LOCATION THROUGH SUBSTRATE (mm 




Figure 4.7 0 y in epoxy fiberglass substrate 



The position where the largest temperature gradients 
in both the X and Y directions are found are shown in Figure 
4.8. These are the points where maximum heat fluxes are 
conducted through the solder/substrate interface. 
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Figure 4.8 Maximum temperature gradients in epoxy fiberglass 
substrate 

2 . High Thermal Conductivity Substrate Packages 

The two substrates chosen in this range are Alumina k 3 
= 18 W/m°C (k 3 /k 1 = 0.12) and Aluminum Nitride (AlN) k 3 = 230 
W/m'C ( k 3 /k 1 = 1.53). Referring back to Figure 4.2 the non- 
dimensional temperature distribution along <p x = .564 for 

Alumina and AlN graph is approximately a vertical straight 
line. Therefore there are no temperature gradients anywhere 
within the package. The reason for no temperature gradients 
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in the substrate is due to the substrate's high thermal 
conductivity, enabling it to diffuse a large heat flux without 
developing a temperature gradient, as shown by Fourier's Law 
of heat conduction, equation 4.7. Figure 4.3 shows the 
largest shear strains at the solder interface have been 
reduced by 40% over the low thermal conductivity substrates 
due to the increased ability of the high thermal conductivity 
substrates to conduct heat away from the chip thereby lowering 
the overall temperature of the package and decreasing the 
thermal expansion. 

There is however no appreciable temperature difference 
between the Alumina and AIN substrate despite the thermal 
conductivity of the Aluminum Nitride being approximately 13 
times as large as the Alumina thermal conductivity. The 
limiting factor for reduction of temperature of the package 
has become the convective thermal resistance, which must be 
lowered in order to lower the package temperature. The 
conductive thermal resistances of the Alumina and AIN 
substrate and convective thermal resistance follow: 



Alumina R cond 



0 . 002 

18 (0 . 00088) (1) 



0.13 °C/W 



0 . 002 



cond 



AIN R, 



230 (0 . 00088) (1) 



0 . 01 °C/W 



conv 200 ( 0 . 001 ) ( 1 ) 

showing that the thermal resistance to heat transfer by 
convection is 38 and 500 times as large as resistance to heat 
transfer by conduction in the Alumina and AlN respectively. 

D. EFFECT OF CONVECTION HEAT TRANSFER ON SYSTEM BEHAVIOR 

This study focuses on the effect of convection cooling on 
the temperature gradients developed within the low thermal 
conductivity substrates and the overall temperature decrease 
of the chip. The electronic packages with high thermal 
conductivity substrates do not develop any temperature 
gradients and exhibit an overall decrease in package 
temperature. The following parameters are held constant; 
volumetric power output (q'" W/m 3 ) of the chip, ambient 
temperature (T„) of the convective fluid, and the following 
materials, silicon chip, Lead/Tin solder and Epoxy Fiberglass 
substrate . 

The governing equation for convection cooling is given by: 

q" = h (r - TJ (4.12) 

where q" is the heat flux in W/m 2 and T is the surface 
temperature. The inverse relation of system temperature to 
convection coefficient is shown in Figure 4.9. This figure 
illustrates the need for ever increasing cooling to maintain 
electronic packages within a temperature tolerance range. An 
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CHIP TEMPERATURE T(°C) 



increase in convective cooling from 100 W/m 2 K to 300 W/m 2 K 
results in a reduction of chip temperature of 90 °C. The same 
increase in convective cooling from 300 W/m 2 K to 600 W/m 2 K 
produces a decrease in chip temperature of only 26 ‘C. 




CONVECTION COEFFICIENT h (W/m*®C) 

Figure 4.9 Convection coefficient vs. chip temperature 

The temperature gradients in both X and Y directions that 
are developed within the low thermal conductivity substrates 
are diminished by increasing the convective heat transfer as 
shown in Figures 4.10 and 4.11. 
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CONVECTION COEFFICIENT h ( W/ m 2 o C) 




TEMPERATURE GRADIENT (°C/mm) 

Figure 4.10 © x in epoxy fiberglass substrate vs. convection 
coefficient 

In the X direction the maximum temperature gradient is 
decreased by 60% by increasing the convective cooling by 250%. 
In the Y direction the maximum temperature gradient is 
decreased by 52% by increasing the convective cooling by 250%. 
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y LOCATION THROUGH SUBSTRATE (mm 




Figure 4.11 6 y in epoxy fiberglass substrate vs. convection 
coefficient 



E. EFFECT OF GEOMETRY ON SYSTEM BEHAVIOR 

The dimensions of the chip and substrate are varied by 
assuming the chip is a percentage of a fixed overall height 
and the solder and substrate make up the remaining percentage. 
We define the non-dimensional term /3 as the ratio of chip 
height to overall package height (H) as in equation 4.13, 
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( 4 . 13 ) 



ft = ^£hip 

H H 

Substrates are approximately kept at a minimum of one half 
a millimeter [Ref. 1]. The following material combination was 
selected; a Silicon chip, a Lead/Tin solder and a Epoxy 
Fiberglass substrate. The following parameters were held 
constant; convective cooling, ambient temperature, the height 
of the air/solder interface and the overall height of the 
package . 

First, a study with a constant heat flux (q") was 
conducted by varying the volumetric heat generation (q'") to 
obtain a flux of 4 x 10 4 W/m 2 . Next the volumetric heat 
generation (q'") was constant which resulted in different heat 
fluxes from the chip. 

1. Constant Heat Flux (q") . 

The result of maintaining a constant heat flux is a 
chip that is a highly concentrated heat source when (3=0.18 and 
a greatly distributed heat source when /3=0.81. In Figure 4.12 
the non-dimensional temperature profile is shown for /3= 0.18 
to 0.81 at a value of <p x = .564. 
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Figure 4 . 12 Non-dimensional temperature for various size chips 
constant heat flux 



When the chip height results in a value of /3=0.18, the 
chip is acting as a concentrated heat source with a volumetric 
heat generation q = 50 MW/m 3 and the temperature of the chip 
is greatest. The larger the chip is to the rest of the 
package the lower the temperature of the overall package is. 
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A partial explanation for this is the amount of surface area 
exposed to convective cooling is smaller when /3 = 0.18 then 
when /3 = 0.81. 

The thermal strain profile is shown if Figure 4.13 for 
the same value of </> x = .564. The thermal strain profiles can 
be explained in much the same way the temperature profiles 
were. That is, when the chip is smallest the largest 
expansion of the package occurs. 




Figure 4.13 Thermal strain profile for various chip heights 
constant heat flux 
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The temperature gradient 0 y is shown in Figure 4.14 
for p values of 0.18 through 0.81. Figure 4.14 shows that 0 y 




TEMPERATURE GRADIENT 0 y (°C/mm) 



Figure 4.14 9 y for various chip heights at constant heat flux 

increases with an increase in p, from 0 y = 36 'C/mm at p = .18 
to 0 y = 57 *C/mm at p = .81. When p = .18 the volumetric heat 
generation q"' is largest at 50 MW/m 3 and provides the least 
surface area exposed to convection. The thermal resistance 
for the substrate is attributable to the temperature gradient 
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developed in the substrate. At /3 = .18 the thermal resistance 



^cond ( P = • 1 8 ) 



3.2 

.16 (.88)1 



22.7 



is large, thereby not allowing much heat flux to conduct 
through the substrate, thereby keeping 0 y small. 

At /3 = .81 the thermal resistance 

i?cond(P = ‘ 81) = .16 ( .88) 1 * 2 ‘ 8 

is very small, allowing for a much larger heat flux to pass 
through the substrate thereby increasing 9 y . Figure 4.15 
shows the temperature gradient 0 x at varying 




Figure 4.15 6 X for various chip heights at constant heat flux 
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positions of <p x for (3 of 0.18 through 0.81. These gradients 
are taken at the top of the substrate corresponding to a value 
of <py = .47. At <p x = .64 0 X is approximately constant at 97 
“C/nun for all values of /3. This can be attributed to the fact 
that no matter what thermal energy is generated by the chip it 
will tend to always flow through the solder and towards the 
upper right convective portion of the substrate. 

In the section of substrate between chips 
corresponding to <p x = .75 and <p x = .84, 0 X increases slightly 
from 31 "C/inin to 52 'C/mm and 14 °C/mm to 27 °C/mm 
respectively. This can be accounted for by looking at the 
substrate equivalent thermal resistances. When /B = .81 the 
resistance is small in the Y direction as compared to the X 
direction as shown by, 



( R cond ^ x 



1.5 

. 16 ( . 4) 1 



23.4 °C/W 



( ^ cond ^ y 



.4 

.16 (.5) 1 



5 °C/W 



Therefore heat transfer would rather take place in the Y 
direction instead of the X direction. When (3 = .18, the 
thermal resistances are large in the Y direction and small in 
the X direction that is, 
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R 



condx 



1.5 

. 16 (3 . 2) 1 



3 °C/W 



R 



condy 



3 . 2 

.16 ( .5) 1 



40 ° C/W 



This change in thermal resistance now favors heat transfer in 
the X direction, which accounts for the increase in 0 x for 
smaller values of /3. 

2. Constant Volumetric Heat Generation (q'") . 

This study shows the effect of using larger power 
chips in increasing sizes. The heat generation of each chip 
is proportional to chip volume, that is, heat generation 
increases with an increase in /3. Figure 4.16 is the 




Figure 4.16 Non-dimensional temperature for various size chips 
constant volumetric power generation 
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normalized temperature profile for (3 = 0.18 to /3 = 0.81 at a 
position of </> x = .564. 

The non dimensional temperature of the package follows 
the identical relationship as did volumetric heat generation 
to (3. As (3 increases so does the non-dimensional temperature 
of the package. If the ratio of heat flux to non-dimensional 
temperature is examined, the increase in non-dimensional 
temperature is found to be non-linear. The following ratios 
q"/T n are given for various values of (3, 



p 


. 18 


.31 


.60 


.73 


.81 


. y n 

— ( v/fliM 


10 


11.2 


16.1 


16.0 


16.0 



The increase in heat flux q”(W/m 2 ) from (3 = .18 to (3 = .31 
is 12 KW/m 2 and from (3 = .6 to (3 = .81 the increase is 14 
KW/m 2 . Even though the increase in q" is 2 KW/m 2 more for the 
later (3's the ratio of q "/T n remains constant. The factors 
that contribute to this are, an increase in the convection 
surface area of the chip and a decrease in the conduction 
thermal resistance of the substrate to allow for a greater 
heat transfer through the substrate. 

Figure 4.17 shows the thermal strain profile at a value of 
</> x = .564. The chip does not have a large increase in thermal 
strain over the range of /3's while the solder has a very large 
increase in thermal strain. This goes back to Table 4.1, 
where the coefficient of thermal expansion for the Pb-Sn 



49 



solder is roughly eleven times greater than that of Silicon. 
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Figure 4.17 Thermal strain profile for various chip heights 
constant volumetric power 



Figure 4.18 shows the variation of the temperature 
gradient Q x at the top of the substrate at various values of 
<p x . At values of <p x = 0.64 and <p x = 0.43 the temperature 
gradient is largest and this corresponds to either side of the 
solder joint. At other values of <p x the temperature gradient 
is diminished greatly. At a value of /3 = 0.81 the temperature 
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gradients are large. They are almost twice as large as 
encountered in any other analysis. For this large a chip 
(large q ") the convection cooling must be increased 
significantly to reduce this temperature gradient. For small 
chip sizes the temperature gradient is very small for any 
value of <p x . 




Figure 4.18 9 X for various chip heights at constant volumetric 
power 
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Figure 4.19 shows the variation in 0 y at the top of 
the substrate where the largest value of the temperature 
gradient within the depth of the substrate occurs. 

When (3 = 0.81 the chip has the greatest heat flux, 
q"=72 KW/m 2 with a correspondingly large 0 y . When /3 = 0.18 the 
chip has the lowest heat flux, q'- 16 KW/m 2 and © y is all but 
non-existent . 




TEMPERATURE GRADIENT 9 (•C/mm) 



Figure 4.19 0 y for various chip heights at constant volumetric 
power 
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3 . Overview Of Geometric Study 

For both the constant heat flux q" (W/m 2 ) and the 
constant volumetric heat generation q'" (W/m 3 ), the temperature 
gradient in the Y direction at the top of the substrate is 
greatest when the chip is largest (/3 = .81). The reasons as 
given previously are the interaction between the conductive 
thermal resistance, the convective thermal resistance, and the 
ratio of volumetric heat generation to available convection 
area . 

If the convective thermal resistance is small as 
compared to the conductive thermal resistance then the 
convective mode will dominate the heat transfer process as the 
heat transfer process is governed by the path of least 
resistance . 

F. EFFECT OF CHIP ENCAPSULATION ON SYSTEM TEMPERATURE 

Some electronic packages are protected from moisture and 
the environment by a protective "encapsulating" coating. For 
this study the encapsulation material EME-1100-T. [Ref . 1] is 
0.1875 mm thick. The effect of encapsulation on the chip 
temperature and gradients developed within the package is 
studied . 

The following parameters will be held constant, volumetric 
heat generation (q'") , convection cooling (h), and ambient 
temperature (T„). These values were the same as in section B 
for comparison purposes. The materials are a Silicon chip, a 
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Lead/Tin Solder, with all four substrates, and the 
encapsulation material. The thermal conductivity of the 
encapsulation material is k = 0.67 W/rn’C. Within the other 
encapsulation materials the thermal conductivity (k) varied 
from 0.67 to 1.97 W/m’C. 

Figure 4.20 is the non-dimensional temperature profile 
taken at a value of <p x = 0.564. The entire package 
experiences less then a 5% temperature change, when low 
thermal conductivity substrates are used and less then 2% 
temperature change when high thermal conductivity substrates 
are used. The encapsulation material causes the surface 
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Figure 4.20 Encapsulated chip non-dimensional temperature 
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NON DIMENSIONAL Y LOCATION <f) 



Even though the encapsulation material has a low thermal 
conductivity its overall thermal resistance is low because its 
less then two tenths of a millimeter thick. 
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cond 
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.67 (2 . 5) 1 



. 11 °C/W 



This is a very low thermal resistance accounting for the small 
change in temperature profile. Figure 4.21 shows the thermal 
strain profile for <p x = .564. 

1.0 

>» 

.8 



.6 



. 4 



.2 



0 

0 .5 1 1.5 2 2.5 

THERMAL STRAIN 6 MT-T )a(x1C>’m/m) 

Figure 4.21 Encapsulated chip thermal strain profile 
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The thermal strain profile has the same negligible 
increases as the normalized temperature profile had except the 
encapsulation material has a very large coefficient of thermal 
expansion (a). 

In summary the effect of encapsulation on system 
temperature and temperature gradients are negligible. 
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V. STRESS ANALYSIS 



A. OVERVIEW OF STRESS ANALYSIS 

The stress analysis conducted is based on using different 
materials in the package and keeping all other factors 
constant. A silicon chip material was used for all cases, 
with two solders Pb-Sn and Au-Sn , and the four substrates 
Epoxy Fiberglass, Polyimide Fiberglass, Alumina and Aluminum 
Nitride. These eight sets of material combinations correspond 
to the thermal strain graphs, shown in Figures 4.3 and 4.4 in 
Chapter IV. 

Figure 4.3 corresponds to the material combinations of a 
silicon chip, the Pb-Sn solder, and the four substrates listed 
above. Figure 4.4 corresponds to the material combinations of 
a silicon chip, the Au-Sn solder, and the four substrates 
listed above. 

The following other factors remain constant, q'"= 20 MW/m 3 
and h = 250 W/m 2 °C. All stresses given with the exception of 
the bending stresses at the centerline of the package are 
along the solder interfaces. The following stresses are 
calculated. First, the bending stresses (a b ), along the 
package centerline. Three stresses are calculated along the 
solder-chip interface, the upper shear stress (t u ), an upper 
bending stress (cr b/U ), and an upper normal stress (a 0 u ) . Three 
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stresses are calculated along the solder-substrate interface, 
the lower shear stress (t L ), a lower bending stress (o bl ), and 
a lower normal stress (a oL ). 

Table 5.1 gives mechanical properties of electronic 
packaging materials, elastic modulus given in GPa. These 
properties are used in the calculation of stresses in 
accordance with the FEM development of chapter three. 

TABLE 5.1 



Material 


Elastic Modulus 


Poisson's Ratio 


Silicon 


13.03 


0.4 


Pb-Sn 


7.4 


7.4 


Pb-Sn 


59.2 


0.4 


EP Fbgls 


35.0 


0.3 


Polyimide 


40.0 


0.3 


Alumina 


262.0 


0.4 


AIN 


339.0 


0.3 



B. STRESS ANALYSIS WITH (Pb-Sn) SOLDER 
1. Bending Stress Along Centerline, o b 

The bending stress at the centerline of the package is 
the lowest of all stresses presented and is shown in Figure 
5.1. The bending stress a b , through the chip is of the order 
of 1 MPa for all four substrates. Through the chip the 
bending stress is tensile for the low thermal conductivity 
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substrates (Epoxy and Polyimide Fiberglasses) and compressive 
for the high thermal conductivity substrates (Alumina and 
AlN) . Through the substrates the bending stresses in the low 
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BENDING STRESS a MPc) 



Figure 5.1 Centerline bending stresses a b 



thermal conductivity substrates is a maximum in compression at 
the top of the substrate ( <p y = .47) at -8 MPa and is a maximum 
in tension at the bottom of the substrate (# y = 0) at 6 MPa. 
The high thermal conductivity substrates have just the 
opposite results, with a maximum in tension at the top of the 
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substrate (<p y = .47) at 2.5 Mpa and maximum in compression at 
the bottom of the substrate (<p y = 0) at -2.5 MPa. Thus the 
bending stresses in the low thermal conductivity substrates 
are about twice as large as the bending stresses in the high 
thermal conductivity substrates . The advantage with respect 
to bending stresses along the centerline of the chip is with 
the high thermal conductivity substrates (alumina and AIN) . 

2. Upper Normal Stress, a Q u 

The upper normal stress distribution along the chip- 
solder interface for all the substrates are compressive with 
a minimum occurring at 6 = .05 and maximum occurring at 6 * 
0.8 as shown in Figure 5.2. For the low thermal conductivity 
substrates the maximum upper normal stress for Epoxy 
Fiberglass is -49 MPa, and the maximum for Polyimide 
Fiberglass is -40 Mpa. Both curves show a decrease in the 
upper normal stress at the edges of the solder due to the fact 
that the edges at 6 = 0.0 and 6=1.0 are free surfaces. For 
the Pb-Sn group of materials the normal stress is the maximum 
stress examined for the low thermal conductivity substrates. 

The upper normal stresses for the high thermal 
conductivity substrates are about -16 MPa and are within 5% of 
each other, decreasing in magnitude slightly at the edges of 
the solder. This shows that the upper normal stresses in the 
high thermal conductivity substrates are about a third of the 
upper normal stress in the low thermal conductivity 
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substrates. We note that the normal stresses are fairly 
uniform from 6 = .25 to 6 = .9. 

In contrast to the chip-solder interface, the normal 
stresses at the solder-substrate interface (i.e. the lower 
normal stresses) are very small on the order of KPa and in 
addition are tensile stresses for all substrate combinations. 
With respect to a 0 u the advantage is again seen with the high 
thermal conductivity substrates (alumina and Ain). 
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Figure 5.2 Upper normal stresses a 
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3. Upper Bending Stress, o b 

The upper bending stress o b u along the chip solder 
interface, is defined as the normal a x stress. However, these 
stresses are presented along the interface, where they are 
greatest, as a function of 6 rather than across the Y 
direction through the solder. These bending stresses are 
compressive for all substrate combinations. For the low 
thermal conductivity substrates the maximum upper bending 
stresses occur at 6 = 0.05 and 0.95 with values of -16 MPa and 
-15 Mpa respectively. For the high thermal conductivity 
substrates maximum upper bending stresses occur at 6 = 0.05 
with -17 Mpa magnitude, and minimums of -8.5 MPa occur at 5 = 



0.8 as shown in Figure 5.3. Again we note that severe 
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gradients occur close to the free edges, and with fairly 
uniform stresses along the center of the chip-solder 
interface. We note these bending stresses fall between the 
bending stresses at the centerline and normal stresses 
reported in the previous sections. With respect to these 
stresses, the advantage (i.e., lower stresses) is obtained for 
the low thermal conductivity substrates (epoxy and polyimide 
fiberglasses ) . 

4. Lower Bending Stress o bl 

The lower bending stresses are compressive in nature 
for all substrate combinations as shown in Figure 5.4. For 
the low thermal conductivity substrates, an immediate decrease 
from about -18 MPa to about -13 MPa, that is about a 30% 
change. Thereafter a gradual decrease to about -11 MPa, that 
is, about 15% change over the remaining domain. The lower 
bending stress for the high thermal conductivity substrates 
remains constant over the entire length of the solder joint at 
-15 MPa. 

In comparison with the upper bending stresses the 
lower bending stresses are greater for the low thermal 
conductivity substrates. The advantage for the lower bending 
stresses goes with the high thermal conductivity substrates 
(alumina and AIN). 
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Figure 5.4 Lower bending stress o b L 
5. Upper Shear Stress r u 

The upper shear stresses at the chip/solder interface 
are positive at 6 = .05 and negative at 6 = .95 for all 
substrate combinations as shown in Figure 5.5. For the low 
thermal conductivity substrates the maximum shear stress is 
negative and occurs at 6 = .95 at -14 MPa and crosses zero 
stress at 5 = .15. For the high thermal conductivity 
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substrates the maximum shear stress is positive and occurs at 
6 = .05 at 11 MPa and crosses zero stress at 6 = .8. The 
upper shear stresses have large gradients at either edge of 
the chip-solder interface, but remain fairly constant through 
the majority of the domain. From 6 = .2 through 6 = .9 the 
upper shear stresses are less than 5 MPa in magnitude. 

These stresses fall in between the upper bending 
stresses and the bending stresses at the centerline. With 
respect to these stresses the advantage goes to the high 
thermal conductivity substrates (alumina and AIN) . 




Figure 5.5 Upper shear stress r u 
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6. Lower Shear Stress r t 

The lower shear stresses are opposite in sign from the 
upper shear stresses, in that the lower shear stresses at the 
solder/substrate interface are negative at 6 = .05 and are 
positive at 6 = .95 for all substrate combinations as shown 
in Figure 5.6. For the low thermal conductivity substrates 
the maximum shear stresses occur at 6 = .05 at about -16 MPa. 
This is comparable in magnitude with the upper shear stresses. 
For the high thermal conductivity substrates the maximum 
positive shear stress occurs at 6 = .95 at 16 MPa while the 
maximum negative shear stress occurs at 6 = .05 at -10 Mpa. 
This stress at 6 = .05 is twice the magnitude as compared to 
the upper shear stresses. Again we note that steep gradients 
occur at the edges of the solder-substrate interface and the 
low thermal conductivity substrates have a gradual decrease 
over the center of the domain. The r L stresses of the low 
thermal conductivity substrates cross zero at 6 = .9. The r L 
stresses for the high thermal conductivity substrates vary 
linearly with 6 between 6 = .2 and .9. These r t stresses 
range from -5 MPa at 6 = .15 to 10 Mpa at 6 = .9 crossing zero 
stress at 6 = .45. We note these lower shear stresses are 
comparable in magnitude to the upper shear stresses. With 
respect to these stresses, there is no clear advantage for 
either set of substrates. 
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Figure 5.6 Lower shear stress r L 



C. STRESS ANALYSIS WITH (Au-Sn) SOLDER 

1. Bending Stress Along Centerline , o b 

The bending stress at the centerline of the package 
remains the lowest of all stresses presented as shown in 
Figure 5.7. The bending stresses at the centerline ( <p x = 0.0) 
for the packages with Au-Sn solder are identical to those with 
Pb-Sn solder. The centerline bending stress is independent of 
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which solder is used. With respect to these stresses, the 
advantage is obtained by the high thermal conductivity 
substrates (alumina and AlN) . 
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Figure 5.7 Centerline bending stress a b 
2. Upper Normal Stress, a Q u 

The upper normal stress distribution along the chip- 
solder interface for all the substrates are compressive in 
nature as shown in Figure 5.8. For the low thermal 
conductivity substrates the maximum for Epoxy Fiberglass is 
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-4 6 MPa and the maximum for Polyimide Fiberglass is -40 Mpa 



and they are constant over a range of .1 < 6 < .9. Both upper 




Figure 5 . 8 Upper normal stress a Q 



normal shear stress distributions decrease at the edges of the 
solder due to the fact that the stress at a free surface is 
zero. The normal stress is once again the maximum stress 
examined for the low thermal conductivity substrates. 

The high thermal conductivity substrates have 
increased by about 20% over the Pb-Sn material package. There 
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are sharp increases at 6 = . 1 and a slight increase to a 
maximum of about -20 MPa at 6 = .9. 

Again we note that in contrast to the chip-solder 
interface, the normal stresses at the solder substrate 
interface are very small on the order of KPa and in addition 
are tensile stresses for all substrate combinations. With 
respect to these stresses, the advantage is obtained for the 
high thermal conductivity substrates (alumina and AIN) and 
with respect to solders, the advantage is obtained for the Pb- 
Sn solder. 

3. Upper Bending Stress, a b u 

The upper bending stress o b along the chip-solder 
interface is defined as the normal o x stress as previously 
mentioned. The upper bending stresses are compressive for all 
substrate combinations as shown in Figure 5.9. For the low 
thermal conductivity substrates maximums occur at 6 = 0.05 and 
0.95 with values of -19 MPa and -22 Mpa respectively. For the 
high thermal conductivity substrates a maximum occurs at 6 = 
0.05 at -39 Mpa. This upper bending stress is twice as large 
for the high thermal conductivity substrates with Au-Sn solder 
then Pb-Sn solder. There is a very large gradient at 6 = 0.05 
where the upper bending stress decreases from -39 MPa to -22 
MPa, that is, about a 45% decrease. For the high thermal 
conductivity substrates the upper bending stress continues to 
decrease at a gradual rate across the remaining domain, 6 = .2 
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to 6 = .9. With respect to these stresses, the advantage is 
obtained for the low thermal conductivity substrates (epoxy 
and polyimide fiberglasses) and with respect to solders, the 
advantage is obtained for the Pb-Sn solder. 
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Figure 5.9 Upper Bending Stress a b u 
4. Lower Bending Stress, o b L 

The lower bending stresses are compressive for all 
substrate combinations as shown in Figure 5.10. For the low 
thermal conductivity substrates a maximum occurs at 6 = 0.05 
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with a value of -39 MPa and a minimum at -7 MPa at 6=0.9. 
This is an increase of 200% over the same package with Pb-Sn 
solder. The lower bending stress for the high thermal 
conductivity substrates remain fairly constant over a range of 
.15 < 6 < .8 with a value of -38 MPa and is a maximum at 6 = 
.95 at -49 MPa. This is an increase of about 260% over the 
same package with Pb-Sn solder. The lower bending stresses 
are the highest stresses presented for the high thermal 
conductivity substrates. With respect to these stresses, the 
advantage is obtained for the low thermal conductivity 
substrates (epoxy and polyimide fiberglasses) and with respect 
to solders, the advantage is obtained for the Pb-Sn solder. 




Figure 5.10 Lower bending stress a b t 
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5. Upper Shear Stress r u 

The upper shear stresses at the chip/solder interface 
are positive at 6 = .05 and negative at 6 = .95 for all 
substrate combinations as shown in Figure 5.11. For the low 
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Figure 5.11 Upper shear stress r u 

thermal conductivity substrates the maximum shear stress is 
negative and occurs at 6 = .95 with a value of -15 MPa and 
crosses zero stress at 6 ~ .15. For the high thermal 
conductivity substrates the maximum shear stress is positive 
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and occurs at 6 



.05 with a value of 16 MPa and crosses zero 



stress at 6 = .8. The maximum upper shear stresses for both 
types of substrates show only a slight increase over the Pb-Sn 
solder packages. With respect to these stresses, there is no 
clear cut advantage for either set of substrates, nor either 
solder . 

6. Lower Shear Stress 

The lower shear stresses are opposite in sign from the 
upper shear stresses, in that the lower shear stresses at the 
solder/substrate interface are negative at 6 = .05 and are 
positive at 6 = .95 for all substrate combinations as shown in 
Figure 5.12. For the low thermal conductivity substrates the 
maximum shear stresses occur at 6 = .05 with a value of -21 
MPa. This represents an increase of 130% over the same 
package with Pb-Sn solder. The lower shear stress then 
decreases to -10 Mpa at 6 = .1 and then slowly decreases to 
about zero at 6 = .95. For the high thermal conductivity 
substrates the maximum shear stresses occur at 6 = .95 with a 
value of 39 MPa. This represents an increase of 240% over the 
same package with Pb-Sn solder. At 6 = .05 the lower shear 
stress is negative at -32 MPa which represents an increase of 
110% over the same package with Pb-Sn solder. With respect to 
these stresses, the advantage is obtained for the low thermal 
conductivity substrates (epoxy and polyimide fiberglasses) and 
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with respect to solders, the advantage is obtained for the Pb- 
Sn solder. 




NORMALIZED X LOCATION ALONG SOLDER 8 = 

L' 



Figure 5.12 Lower Shear Stress r t 

D. SUMMARY OF STRESS ANALYSIS 

The following table represents the maximum stresses (given 
in MPa) for each material combination studied. Various 
comparisons are possible. 
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TABLE 5.2 MAXIMUM STRESSES 



EP FBGLS 



POLY FBGLS 



ALUMINA 



AIN 





Pb-Sn 


Au-Sn 


a b 


-8 


-8 


a 0 


-49 


-4 6 


a b,u 


-16 


-20 


CT b,l 


-19.5 


-39 


r u 


-14.5 


-16 


r l 


-16 


-21 


a b 


-7.5 


- 7.5 


a o 


-41 


-40 


a b,u 


-14 


-18 


a M 


-16.5 


-34 


r u 


-13 


-14 




-15 


-20 


a b 


-2.5 


-3 


a o 


-16 


-20 


a b,u 


-17 


-38 


CT b,l 


-15 


-45 


r u 


11 


16 


r L 


15 


40 


°b 


-2.5 


-3 


a o 


-15 


-19 


CT b,u 


-17 


-38 


a M 


-14 


-43 


r u 


11 


' 16 


r l 


15 


40 
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For the Epoxy fiberglass substrate packages there is 
little difference in centerline bending a b , normal a 0 , and 
upper shear r u stresses and only a slightly larger upper 
bending a b u , and lower shear stresses for Au-Sn versus Pb- 
Sn solder. For this substrate the lower bending a bl stresses 
are twice as great for Au-Sn versus Pb-Sn solder. The normal 
stresses are the largest stresses for this substrate for both 
Pb-Sn and Au-Sn solders. The most significant stresses with 
regard to deforming the solder joint are the lower bending and 
lower shear. The upper normal stress is not significant 
because it tends to compress a solder crack rather then 
propagate it. 

For the Polyimide fiberglass substrate packages there is 
little difference in centerline bending a b , normal a 0 , and 
upper shear r u stresses and only a slightly larger upper 
bending a b u , and lower shear t x stresses for Au-Sn versus Pb- 
Sn solder. For this substrate the lower bending a b x stresses 
are twice as great for Au-Sn versus Pb-Sn solder. The normal 
stresses are the largest stresses for this substrate for both 
Pb-Sn and Au-Sn solders. The most significant stresses with 
regard to deforming the solder joint are the lower bending and 
lower shear. The upper normal stress is not significant 
because it tends to compress a solder crack rather then 
propogate it. The Polyimide fiberglass packages have stress 
distributions similiar to Epoxy fiberglass, but tend to be 
slightly lower in magnitude. 
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For the Alumina substrate only the centerline bending o b 
stresses are similiar for both solders. The upper normal a Q 
and upper shear r u stresses are slightly larger for Au-Sn 
versus Pb-Sn solders. The upper bending o b u stress is 225% 
larger, lower bending a b L stress is 300% larger, and the lower 
shear r t stress is 265% larger for the Au-Sn versus Pb-Sn 
solder. The upper bending stress is largest for the Pb-Sn 
solder while the lower bending stress is largest for the Au-Sn 
solder. The most significant stresses with regard to 
deforming the solder joint are the upper bending, lower 
bending and lower shear. 

For the AIN substrate only the centerline bending a b 
stresses are similiar for both solders. The upper normal a 0 
and upper shear r u stresses are slightly larger for Au-Sn 
versus Pb-Sn solders. The upper bending a b u stress is 225% 
larger, lower bending a b L stress is 305% larger, and the lower 
shear r t stress is 265% larger for the Au-Sn versus Pb-Sn 
solder. The upper bending stress is largest for the Pb-Sn 
solder while the lower bending stress is largest for the Au-Sn 
solder. The most significant stresses with regard to 
deforming the solder joint are the upper bending, lower 
bending and lower shear. The AlN packages have stress 
distributions which are approximately less than 5% different 
from Alumina. 

Based on the stress results from this investigation the 
following ranking of material combinations are presented based 
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on obtaining the lowest stresses in the package, 
chip is used in all packages. 



The silicon 



1. 


Pb-Sn 


with 


Alumina or Ain 


2. 


Pb-Sn 


with 


Polyimide fiberglass 


3. 


Pb-Sn 


with 


Epoxy fiberglass 


4. 


Au-Sn 


with 


Polyimide fiberglass 


5. 


Au-Sn 


with 


Epoxy fiberglass 


6. 


Au-Sn 


with 


Alumina or Ain 


Bending 


K) 


and Shear (t) stresses at 



substrate interface can be significant. The lower inside 
corner of the solder-substrate interface (<S = 0.0) is the 
point where failure is most likely to occur due to bending and 
shear stresses. 

This investigator recommends the following topics for 
further research. Alter the thermal code for transient 
temperature response of the package. The result of this would 
be to track stress in time. In this way the changes in 
temperature and stress can be tracked in time after the 
electronic device has been turned on. Conduct a parametric 
study of a general tri-material package varying the materials 
coefficient of thermal expansion and modulus of elasticity to 
better understand the effect on stresses produced in the tri- 
material package. 
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APPENDIX A 



PROGRAM HEAT 

********************************************************************** 
* * 

* This F.E.M. program solves for the Temperature distribution in * 

* electronic packages. The governing differential equation is * 

* POISSON's equation. More specifically: * 

* * 

* * 



When symmetry is 
are used at these 



envoked homogeneous Neumann boundary 
locations. The other boundaries are 



conditions 
Cauchy . 



Listing of variables: 

= local nodal point 



aai j 
ALPHA 
AREA 
bi 

bbi j 

BIGA 

BIGF 

BDRYVAL 

BLOW 

cauchf i 

ci 

DELI , 2 , 3 

DELMAX 

DIVSUB 

DIVX 

DIVY 

DIVXAIR 

DIVYAIR 

ffl 

GEN 

HLOW 

HUP 

HVERT 

I MAT 

I COR 

KEL 

KINVS 

NBCLOW 

NBCUP 

NBCVERT 

NBDLEFT 

NBDRIT 

NBDUP 

NBDLOW 

NCOL 

NROW 

NSUBEL 

NAIREL 



values of BIGA matrix 



vector of element areas 

triangular coordinate property of each element 
contribution from Cauchy boundary to BIGA matrix 
global left hand matrix 
right hand side vector 
value of boundary point 

( BUP , BVERT ) inverse of the absolute temp of the fluid 
contribution from Cauchy boundary to BIGF vector 
triangular coordinate property of each element 



surface temps 



grid 

grid 



difference between assumed and calc 
maximum of dell or del2 or del3 

# of rows of substrate 

# of divisions in x direction for a 

# of divisions in y direction for a 

# of columns of air in grid 

# of rows of air(solder) in grid 
contribution of source term to BIGF vector 
if equal to 1 then internal heat generation 
lower side heat transfer coefficient 
upper side heat transfer coefficient 

right side heat transfer coefficient 
identifies what material properties go with 
establishes correspondence between local and 
therm, conductivity of element 
kinematic viscosity of the fluid 



density 

density 



is involved 



* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

★ 

* 

* 

* 

★ 



which element* 
global coor.* 



lower boundary 
upper boundary 
right boundary 



is a Cauchy 

is a Cauchy 

is a Cauchy 

left boundary is a Direchlet 
right boundary is a Direchlet 

upper boundary is a Direchlet 

lower boundary is a Direchlet 

number of columns for grid subroutine 
number of rows for grid subroutine 

# of substrate elements 

# of columns of air in grid 



if 

if 

if 

if 

if 

if 

if 



then 

then 

then 

then 

then 

then 

then 



* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 
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&RSOEL 

DOT 

MjLOW 

\LUP 

\LVERT 

-10 

30EL 

PHT 

PHTEL 

\MB 

SMP 

FLLOW 

FLUP 

FLVERT 

SLOW 

SUP 

3VERT 

IERCON 

POS 

POS 

******* 



* # of air & solder elements 

- volumetric heat generation of electronic device 

- Raleigh # of the fluid at the lower surface 

■» Raleigh # of the fluid at the upper surface 

=■ Raleigh # of the fluid at the vertical surface 

* density of 4 different materials 

= density of element 

= specific heat of 4 different materials 
= specific heat of element 

- ambient temperature 

= solution vector of temperatures 
= fluid temperature at lower surface 
= fluid temperature at upper surface 
■ fluid temperature at vertical surface 

- average surface temperature of lower surface 
= average surface temperature of upper surface 

- average surface temperature of vertical surface 
= thermal conductivity of 4 different materials 

- x position of system nodal point 
= x position of new grid density 

= y position of system nodal point 

- y position of new grid density 
***************************************************** 



* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

******* 



INCLUDE ' HEATCOMN. FOR' 

***** ensure the proper amount of workspace is allocated *** 
***** on problems with a large quantity of unknowns a ****** 
***** FORTRAN STOP may be encountered. the amount of ****** 
***** work space allocated in RWKSP and IWKIN must be ****** 
***** increased 

COMMON/WORKSP/RWKSP 
REAL RWKSP (900000) 

CALL IWKIN (900000) 



****** ensure the proper data files are being read ********** 

OPEN ( 7 , FI LE= ' HEAT . OUT ' , STATUS= ' NEW ' ) 

OPEN ( 8 , FI LE= ' TESTl . OUT ' , STATUS= ' NEW ' ) 

OPEN ( 9 , FI LE= ' TEST2 . OUT ' , STATUS- ' NEW ' ) 

reads given data and creates mesh 
CALL INPUT2 
CALL GRIDPLUS 

1 creates biga matrix and bigf vector 
CALL HEATMAT 

replaces Galerkin equation with direchlet bdry temperatures 
1 modifies equations so that heat flux at boundaries of material 
1 interfaces are equal 
CALL HEATMOD 

1 solves simultaneous system of equations 
CALL LSARG (NSNP, BIGA, NSNP, BIGF, 1 , TEMP) 

creates output file 
CALL OUTPUT 



END 
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SUBROUTINE INPUT2 



* * * * 
★ * ★ ★ 



★ ★ ★ ★ 
* * * * 
* ★ ★ ★ 
★ * * * 



C 



* * * ★ 
* * * * 



************************************************************** 
Reads in data to construct thermal grid 



INCLUDE 'HEATCOMN. FOR' 

OPEN( 11, FILE- 'da ta8.DAT' , STATUS- ' OLD ' ) 
OPEN ( 21 , FILE-' TEST . OUT ' , STATUS- ' NEW ' ) 



READ( 11 , * ) 
READ( 11 , * ) 
READ ( 11 , * ) 
READ ( 11 , * ) 
READ ( 11 , * ) 
READ ( 11 , * ) 
READ( 11 , * ) 
READ( 11 , * ) 
READ( 11 , * ) 
READ( 11 , * ) 
READ( 11 , * ) 
READ( 11 , * ) 



NCOLl ,NCOL2 ,NROW 
( XPOS ( I ) , 1=1 ,NDIFX) 

( YPOS ( I ) , 1=1 ,NDIFY) 

( XINCOL ( I ) , 1=1 , NDI FX-1 ) 

( YINCRROW ( I ) , 1=1 ,NDIFY-1 ) 

( THERCON ( I ) , RHO ( I ) , SPHT( I ) , 1=1 , NMATL ) 
DIVSUB , DIVXAIR , DIVYAIR , DIVXSOL 
NBCLOW , NBCUP , NBCVERT , NBCSUB , TAMB 
HLOW, HUP , HVERT , HSUB 
GEN , QTPR 

NBDUP , NBDRDEV , NBDUSUB , NBDRIT , NBDLOW 
TUP , TRTDEV , TUPSUB , TRIGHT , TLOW 



END 

SUBROUTINE GRIDPLUS 



********************************************************** 
overlays a grid on a square or rectangle with up to nine 
different mesh densities in X and Y directions , sets up 
geometric and material properties for the Galerkin eq. 

INCLUDE 'HEATCOMN. FOR' 

open( 21, file-' test. out' , status- 'new' ) 

DATA KEL , XVALUE , YVALUE /NEL*0,0,0/ 

RBMAXl = NCOLl + 1 
RBMAX2 = NCOL2 + 1 
RTMAX = NROW + 1 
XCOUNT = 1 
NN =1 

RBMAX = RBMAXl 



*********************************************************** 
assigns X and Y coordinates to system Global nodal points 

DO 100 I = 1 , NSNP 

X(I) = XVALUE 
Y ( I ) = YVALUE 



I F( XVALUE.LT. XPOS ( 2 ) . AND . XVALUE . GE . XPOS ( 1 ) ) THEN 



XINCCOL = XINCOL ( 1 ) 
ELSE I F( XVALUE.LT. XPOS ( 3 
XINCCOL = XINCOL ( 2 ) 
ELSEI F( XVALUE.LT. XPOS ( 4 
XINCCOL = XINCOL ( 3 ) 
ELSEI F( XVALUE.LT. XPOS ( 5 
XINCCOL - XINCOL ( 4 ) 
ELSEI F( XVALUE.LT. XPOS (6 
XINCCOL - XINCOL ( 5 ) 



.AND. XVALUE. GE. XPOS ( 2 ) ) THEN 
. AND . XVALUE . GE . XPOS ( 3 ) ) THEN 
.AND. XVALUE. GE. XPOS ( 4 ) ) THEN 
.AND. XVALUE. GE.XPOS( 5) ) THEN 



* * * * 



82 



ELSEI F ( X VALUE . LT . XPOS ( 9 ) . AND . XVALUE . GE . XPOS ( 8 ) ) THEN 
XINCCOL = XINCOL ( 8 ) 

ELSE 

XINCCOL = XINCOL ( 9 ) 

ENDIF 

XVALUE = XVALUE + XINCCOL 

I F ( YVALUE . LT . YPOS ( 2 ) . AND . YVALUE . GE . YPOS ( 1 ) ) THEN 
YINCROW = YINCRROW ( 1 ) 

ELSEI F( YVALUE.LT. YPOS ( 3 ) . AND . YVALUE . GE . YPOS ( 2 ) ) THEN 
YINCROW = YINCRROW ( 2 ) 

ELSEI F( YVALUE . LT . YPOS ( 4 ) . AND . YVALUE . GE . YPOS ( 3 ) ) THEN 
YINCROW - YINCRROW ( 3 ) 

ELSEI F( YVALUE.LT. YPOS ( 5) . AND . YVALUE . GE . YPOS ( 4 ) ) THEN 
YINCROW = YINCRROW ( 4 ) 

ELSEIF( YVALUE. LT.YPOS( 6) . AND . YVALUE . GE . YPOS ( 5) ) THEN 
YINCROW = YINCRROW ( 5 ) 

ELSEI F ( YVALUE . LT . YPOS ( 7 ) . AND . YVALUE . GE . YPOS ( 6 ) ) THEN 
YINCROW = YINCRROW ( 6 ) 



ELSEI F ( YVALUE . LT . YPOS ( 8 ) . AND . YVALUE . GE . YPOS ( 7 ) ) THEN 
YINCROW - YINCRROW ( 7 ) 

ELSEI F ( YVALUE . LT . YPOS ( 9 ) . AND . YVALUE . GE . YPOS ( 8 ) ) THEN 
YINCROW = YINCRROW ( 8 ) 

ELSE 

YINCROW = YINCRROW ( 9 ) 

ENDIF 

begins numbering process of assembly after substrate 
I F ( I . GT . ( ( NCOLl+1 ) * ( DIVSUB+1 ) ) ) THEN 
RBMAX = RBMAX2 
ENDIF 

I F ( XCOUNT . GE . RBMAX ) THEN 
XVALUE = 0.0 

YVALUE = YVALUE + YINCROW 
XCOUNT = 0.0 
ENDIF 



XCOUNT = XCOUNT + 1 
CONTINUE 

************************************************************* 
creates correspondence table for Local-Global nodal points 

DO 200 IROW - 1 , DIVSUB + 1 
KK - ( I ROW— 1 ) * ( RBMAXl ) 

DO 210 ICOL - 1 , NCOLl 



I EL - ( (IROW-l)*(2*NCOLl) ) + (2*ICOL)-l 



JEL 



IEL + 1 



I COR ( IEL , 1 ) 
I COR ( IEL, 2 ) 
I COR ( IEL, 3 ) 
ICOR( JEL, 1 ) 
ICOR( JEL, 2 ) 
I COR ( JEL , 3 ) 



ICOL 

ICOL 

ICOL 

ICOL 

ICOL 

ICOL 



KK 

NCOLl + 
NCOLl + 
KK 

KK + 1 
NCOLl + 



KK 

KK 



2 

1 



KK + 2 
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210 CONTINUE 

200 CONTINUE 

**** correspondence starting two rows above substrate 
KK = ( RBMAXl* ( DIVSUB+1 ) ) 

I EL = ( NCOLl*DIVSUB*2 ) + (NCOL2*2) + 1 
DO 220 IROW - l,NROW-( DIVSUB+1) 

DO 230 ICOL = 1 , NCOL2 

JEL = I EL + 1 



I COR ( I EL , 1 ) 


= 


ICOL 


+ 


KK 








I COR ( I EL , 2 ) 


= 


ICOL 


+ 


NCOL2 + 


KK 


+ 


2 


I COR( I EL , 3 ) 


= 


ICOL 


+ 


NCOL2 + 


KK 


+ 


1 


ICOR( JEL, 1 ) 


= 


ICOL 


+ 


KK 








ICOR( JEL, 2 ) 


= 


ICOL 


+ 


KK + 1 








ICOR( JEL, 3 ) 


= 


ICOL 


+ 


NCOL2 + 


KK 


+ 


2 



I EL = I EL + 2 

230 CONTINUE 

KK «* KK + RBMAX2 
220 CONTINUE 

*************************************************************** 
**** assigns material properties to elements 

NSUBEL = NCOLl * DIVSUB * 2 
NAIREL = 2 * DIVXAIR 

NARSOEL = 2 * DIVYAIR * NCOL2 

LCOUNT = 1 

NCOUNT - 1 

MCOUNT = 1 

DO 300 I « 1 , NEL-1 , 2 
I F ( I . LT. NSUBEL) THEN 
**** substrate (single material) **** 

I MAT ( I ) = 1 

I MAT ( 1+1 ) = 1 

**** substrate (multiple material) **** 

ELSEI F ( I . GT . NSUBEL . AND . I . LT . ( NARSOEL+NSUBEL ) ) THEN 
I F ( NCOUNT . LE . DIVXAIR ) THEN 
★ * ★ * sir * * * * 

I MAT ( I ) = 2 

I MAT ( I +1 ) - 2 

ELSEI F ( NCOUNT . LE . DI VXAIR+DIVXSOL ) THEN 
**** solder **** 

I MAT ( I ) = 3 

IMAT( 1+1 ) = 3 
**** coating/cover **** 

ELSE 

I MAT ( I ) = 5 

I MAT ( I +1 ) = 5 
ENDIF 

NCOUNT = NCOUNT + 1 
I F ( NCOUNT . GT . NCOL2 ) THEN 
NCOUNT = 1 
ENDIF 

ELSE I F ( I .GT. NARSOEL+NSUBEL. AND. 

: I . LT . NEL- ( NCOL2 * 2 ) ) THEN 

I F ( MCOUNT . LE . NCOL2-1 ) THEN 
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device **** 

I MAT ( I ) = 4 

IMAT( 1+1 ) = 4 
coating/cover **** 

ELSE 

I MAT ( I ) = 5 

I MAT ( I + 1 ) = 5 
ENDIF 

MCOUNT = MCOUNT + 1 
IF(MCOUNT.GT.NCOL2 )THEN 
MCOUNT = 1 
ENDIF 
ELSE 

coating/cover **** 

I MAT ( I ) = 5 

IMAT( 1 + 1 ) = 5 
ENDIF 
CONTINUE 



DO 310 I = 
II 

KEL ( I ) 
RHOEL ( I 



1 , NEL 

=■ I MAT ( I ) 

= THERCON (II) 
= RHO (II) 



SPHTEL ( I ) = SPHT(II) 

I F ( I .GT. ( NSUBEL+NARSOEL ) ) THEN 
QDOTEL ( I ) = QTPR 
ENDIF 
CONTINUE 



*********************************************************** 
assigns boundary conditions to elements at the convective 
(CAUCHY) boundaries. 

DO 400 I = 2 , 2 *NCOLl , 2 

assigns Cauchy boundary to lower side of rectangle 
I F ( NBCLOW . EQ . 1 ) THEN 
NBDRY(I) = 1 
ENDIF 
CONTINUE 

DO 410 I = ( NEL- ( 2 *NCOL2 ) ) + 1 , NEL-1 , 2 

assigns Cauchy boundary to upper side of rectangle 
I F ( NBCUP . EQ . 1 ) THEN 
NBDRY(I) - 2 
ENDIF 
CONTINUE 



DO 420 I = NSUBEL- ( (NCOLl-NCOL2)*2)+l, NSUBEL-1 , 2 

issigns Cauchy boundary to upper tip of substrate 
I F ( NBCSUB . EQ . 1 ) THEN 
NBDRY(I) « 4 
ENDIF 
* CONTINUE 
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DO 430 I = NSUBEL+ ( NCOL2 * 2 ) ,NEL,NCOL2*2 

*** assigns Cauchy boundary to vertical side of rectangle 
I F ( NBC VERT . EQ . 1 ) THEN 
NBDRY(I) = 3 
ENDIF 
430 CONTINUE 

************************************************************** 

**** assigns Direchlet boundary to system Global nodal points 
* * * * 

**** assigns Direchlet boundary to lower side of substrate 

DO 500 I = 1 , RBMAXl 

IF( NBDLOW. EQ. 1 ) THEN 
NDIRECH(I) = 1 
BDRYVAL ( I ) = TLOW 
ENDIF 

500 CONTINUE 

**** assigns Direchlet boundary to right hand side of substrate 

DO 510 I = RBMAXl *2, RBMAXl *( DI VSUB+1 ), RBMAXl 
I F ( NBDRI T . EQ . 1 ) THEN 
NDIRECH ( I ) = 1 



BDRYVAL ( I ) = TRIGHT 
ENDIF 
510 CONTINUE 

**** assigns Direchlet boundary to upper side of device 

DO 520 I = NSNP , NSNP— NCOL2 , — 1 
I F ( NBDUP . EQ . 1 ) THEN 
NDIRECH(I) = 1 
BDRYVAL ( I ) = TUP 
ENDIF 
520 CONTINUE 

**** assigns Direchlet boundary to right side of device 

III = NSNP- ( RBMAX2 * ( NROW-DI VSUB-1 ) ) 

DO 530 I = NSNP-RBMAX2 , 1 1 1 , -RBMAX2 
I F ( NBDRDEV . EQ . 1 ) THEN 
NDIRECH(I) - 1 
BDRYVAL ( I ) = TRTDEV 
ENDIF 
530 CONTINUE 

**** assigns Direchlet boundary to upper side of substrate 

DO 540 I = ( RBMAXl *DI VSUB ) +RBMAX2 , RBMAXl *(DIVSUB+1) 

I F ( NBDUSUB . EQ . 1 ) THEN 
NDIRECH ( I ) = 1 
BDRYVAL ( I ) = TUPSUB 
ENDIF 
540 CONTINUE 

END 
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SUBROUTINE HEATMAT 



***************************** 
calculates the BIGA matrix a 
system of simultaneous equat 
for the main operator ( Del "2 ) 
forcing function to Poisson' 

INCLUDE ' HEATCOMN. FOR' 



DATA BIGA, BIGF / 

DATA bb21 , bb22 , bb23 , bb33 /0 
DATA cauchf 1 , cauchf 2 , f f 1 /0 

DO 10 I EL = 1 , NEL 



11 = 

12 = 
13 = 

Bl = 
B2 - 
B3 - 

Cl = 
C2 - 
C3 = 

AREA( 

write 

aall 

aal2 

aal 3 

aa21 

aa22 

aa23 

aa31 

aa32 

aa33 



ICOR ( I 
I COR ( I 
I COR ( I 

Y(I2) 

Y(I3) 

Y(I1) 

X(I3) 
X(I1) 
X( 12) 

I EL ) - 

( 8 ,*) 

= KEL ( 
= KEL ( 
= KEL ( 
= KEL ( 
= KEL ( 
= KEL ( 
= KEL ( 
= KEL ( 
= KEL ( 



EL,1) 
EL, 2) 
EL, 3) 

- Y ( I 

- Y ( I 

- Y ( I 

- X(I 

- X(I 

- X(I 

( ABS 
iel , a 
I EL ) * 
IEL) * 
IEL) * 
IEL) * 
IEL) * 
IEL) * 
IEL) * 
IEL) * 
IEL) * 



3) 

1) 

2 ) 

2 ) 

3) 

1) 



(Cl* 
rea ( 
( Bl* 
( Bl* 
(Bl* 
( B2* 
(B2* 
( B2* 
(B3* 
(B3* 
( B3 * 



***************************** 
calculation of Cauchy bounda 
vector 



**************************** 
nd BIGF vector that are the 
ions to be solved; accounts 
, Cauchy boundary, and the 
s equation 



NDIM* 0 , NSNP* 0 / 
. 0 , 0 . 0 , 0 . 0 , 0 . 0 / 
. 0 , 0 . 0 , 0 . 0 / 



B2-C2 *Bl ) )/2 . 
iel ) 

*2 + Cl**2 )/( -4 . *AREA( IEL) ) 

B2 + Cl*C2)/(-4. *AREA( IEL ) ) 
B3 + Cl*C3 )/(-4 . *AREA( IEL) ) 
Bl + C2*Cl)/(-4. *AREA( IEL ) ) 
*2 + C2**2 )/(-4 . *AREA( IEL) ) 

B3 + C2*C3)/(-4. *AREA( IEL) ) 
Bl + C3*Cl)/(-4 . *AREA( IEL) ) 
B2 + C3*C2 )/(-4 . *AREA( IEL) ) 
*2 + C3**2 )/(-4 . *AREA( IEL) ) 



***************************** 
ry to BIGA matrix and BIGF 



IF(NBDRY( IEL) . NE . 0 ) THEN 

boundary with local nodal points 1 - 2 on bottom side 
I F ( NBDRY ( IEL) . EQ . 1 ) THEN 

bb22 = ( H LOW* ABS ( C3 ) )/-3 . 0 
bb21 = ( BLOW* ABS ( C3 ) )/-6 . 0 
cauchf 1 = ( TAMB*HLOW*ABS ( C3 ) )/— 2 . 0 
boundary with local nodal points 2 - 3 on upper side 
ELSEI F ( NBDRY (IEL) . EQ . 2 ) THEN 

bb33 = (HUP*ABS(Cl) ) /— 3 . 0 
bb23 = ( HUP*ABS( Cl ) )/-6 . 0 
cauchf 2 = ( TAM B* HUP* ABS ( Cl ) )/— 2 . 0 
boundary with local nodal points 2 - 3 on right side 
ELSEI F( NBDRY ( IEL) . EQ. 3 )THEN 

bb33 = ( H VERT* ABS ( Bl ) )/-3 . 0 
bb23 = ( HVERT*ABS( Bl ) )/-6 . 0 
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cauchf 2 
ELSE 

bb3 3 
bb2 3 
cauchf 2 
ENDIF 
ENDIF 



( TAMB*HVERT*ABS ( Bl ) )/-2.0 

( HSUB*ABS( Cl ) ) /- 3.0 
( HSUB *ABS ( Cl ) )/-6 .0 
( TAMB*HSUB*ABS ( Cl ) )/-2 .0 



******************* 
**** the contributi 
**** bbll=bb22 bbl2 



B I GA ( I 
BIGA( I 
BIGA( I 
BIGA( I 
BIGA( I 
BIGA( I 
BIGA( I 
BIGA( I 
BIGA( I 



**************************************** 
on of the Cauchy boundary is symmetric 
= bb21 bb22=bb3 3 bb23=bb32 



1,11) 


= 


BIGA( 


1,12) 


= 


BIGA( 


1,13) 


= 


BIGA( 


2,11) 


= 


BIGA( 


2,12) 


= 


BIGA( 


2,13) 


= 


BIGA( 


3,11) 


= 


BIGA( 


3,12) 


= 


BIGA( 


3,13) 


= 


BIGA( 



11,11) 


+ 


aall 


11,12) 


+ 


aal2 


11,13) 


+ 


aal3 


12,11) 


+ 


aa21 


12,12) 


+ 


aa2 2 


12,13) 


+ 


aa23 


13,11) 


+ 


aa31 


13,12) 


+ 


aa32 


13,13) 


+ 


aa 33 



+ bb22 
+ bb21 

+ bb21 

+ bb22 + bb3 3 
+ bb2 3 

+ bb2 3 
+ bb3 3 



bb22 = 0.0 
bb21 = 0.0 
bb2 3 = 0.0 
bb3 3 = 0.0 



************************************************************ 
**** US ed for Poisson's equation (internal heat generation) 

IF( GEN. EQ. 1.0) THEN 

I F ( I MAT ( I EL) .EQ. 4 ) THEN 

ffl = -AREA( IEL) *QTPREL( IEL)/3 .0 
ENDIF 

ENDIF 



*************************************************************** 
**** the two terms that comprise the BIGF matrix are from the 
**** Cauchy boundary condition and the source term of the chip 
**** a lumped approximation is used for both cases; cauchfl is 
**** for a 1-2 boundary edge cauchf2 is for a 2-3 boundary edge 



BIGF ( I 1 ) = BIGF ( 
B IGF ( 1 2 ) = BIGF ( 
BIGF ( 13 ) = BIGF ( 



ID 


+ 


ffl 


+ 


cau 


12) 


+ 


ffl 


+ 


cau 


13) 


+ 


ffl 


+ 


cau 



chf 1 
chf 1 
chf 2 



+ cauchf2 



cauchfl = 0.0 
cauchf2 =* 0.0 
ffl - 0.0 



10 CONTINUE 
END 



SUBROUTINE HEATMOD 



****************************************************************** 
**** modifies BIGA matrix and BIGF vector if a system Global nodal 
**** point is a Direchlet boundary ie. replaces the Galerkin eq. 
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INCLUDE ' HEATCOMN. FOR 

modifies system of equations for Direchlet boundary points 
DO 10 I ■ 1 , NSNP 

IF ( NDIRECH ( I ) .EQ. 1 ) THEN 
BIGA( 1,1) = 1.0 

DO 20 J = 1,1-1 

B I GA ( I , J ) = 0.0 
CONTINUE 

DO 30 JJ = I +1 , NSNP 
B I GA ( I , J J ) = 0.0 
CONTINUE 

BIGF(I) = BDRYVAL ( I ) 

ENDIF 

CONTINUE 

END 

SUBROUTINE OUTPUT 

**************************************************************** 
creates output file of temperature distribution and a data file 
for graphing the temperature contours 

INCLUDE ' HEATCOMN . FOR ' 

OPEN ( 10, FILE-' HEAT. DAT' , STATUS= ' NEW ' ) 
write(7,*) ' Temperatures through solder' 

write ( 7 , * ) '8', temp( 8) , '23' ,temp( 23) , '38' , temp( 38 ) , ' 53 ' , temp( 53 ) 

write ( 7, * ) '68' , temp (68) , '83' , temp (83) , '93' , temp( 93 ) , 

: ' 103 ' , temp( 103 ) 

write ( 7 , * ) '113', temp (113), '123', temp( 123 ) , ' 133 ' , temp( 133 ) 

energy balance 

qconv = 0.0 

do 200 i = 126,134 

qconv-qconv+ ( hup* ( ( ( temp( I ) +temp( 1+1 ) )/2 )-tamb )*(X(I+l)-x(l))) 
continue 

do 210 i - 1,14 

qconv=qconv+ ( hlow* ( ( ( temp( I )+temp( 1+1 ) )/2 )-tamb) *(X(I+l)-x(l))) 
continue 

do 220 i = 70,74 

qconv=qconv+ ( hsub* ( ( ( temp ( I ) +temp ( 1+1 ) )/2 ) -tamb )*(X(I+l)-x(I) ) ) 
continue 

do 230 i = 85,125,10 

qconv=qconv+( hve r t * ( ( ( temp( I )+temp( 1 + 10 ) )/2 )-tamb) * 
s (y( 1+10 )-y( I ) ) ) 

continue 

qconv=qconv+( hver t* ( ( ( temp( 70 ) +temp ( 85 ) )/2 ) -tamb ) * ( y ( 85 )-y ( 70 ) ) ) 

WRITE ( 7 , * ) 

write(7,*) ' Heat loss to convection = ', qconv, 'W/M' 

POWER = QPRT*0 . 0025*0 . 002 

WRITE ( 7 , * ) 'The power output of the device - ' , powe r , 'W/M' 

END 
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APPENDIX B. 



* * * * 



PROGRAM FILTER 

REAL X ( 135) ,Y(135),T(135) , TEMPT (200) , TEMPB(200) 

OPEN ( 12 , FILE= ' bh3 . out ' , STATUS® ' OLD ' ) 

OPEN ( 13, FILE=' bh3f.DAT' , STATUS® ' NEW ' ) 

READ ( 12, * ) (X( I ) , Y( I ) ,T( I ) , 1=1 , 135) 

calculates intermediate temperatures ******* 



TA = 


( T ( 61 ) -T( 4 6 ) ) * . 7 5 + 


T( 46 ) 


TB = 


( T( 61 ) +T( 46 ) )/2.0 




TC = 


( T ( 6 2 ) -T ( 4 7 ) ) * . 7 5 + 


T( 47 ) 


TD = 


(T(66)+T( 47) )/2 . 0 




TE - 


(T( 63 )-T( 48 ) ) * . 75 + 


T ( 4 8 ) 


TF = 


( T( 6 3 ) +T( 4 8 ) )/2 . 0 




TG = 


(T( 64 )-T( 49 ) ) * . 75 + 


T( 49 ) 


TH = 


( T ( 64 )+T( 49 ) )/2 . 0 




TI = 


( T( 65 )-T( 64 ) ) * . 75 + 


T( 64 ) 


TL = 


( T( 50 ) -T( 49 ) ) * . 7 5 + 


T( 49 ) 


TJ = 


( TI-TL ) * . 7 5 + TL 




TK = 


( TI+TL)/2 . 0 




TM = 


( T( 6 5 ) -T ( 50 ) ) * . 7 5 + 


T ( 50 ) 


TN = 


( T ( 6 5 ) +T ( 5 0 ) )/2 . 0 




TO = 


( T ( 6 5 ) +T ( 66 ) )/2 . 0 




TR - 


( T ( 51) -T( 50 ) ) * . 75 + 


T ( 50 ) 


TP = 


( TO— TR ) * . 7 5 + TR 




TQ = 


( TO+TR ) /2 . 0 




TS = 


( T ( 6 5 )+T( 66 ) )/2 . 0 




TV = 


( T ( 50 ) +T ( 51) )/2 . 0 




TT = 


( TS-TV ) * . 7 5 + TV 




TU = 


( TS+TV ) /2 . 0 




TW = 


( T( 66 )-T( 51 ) ) * . 75 + 


T ( 51) 


TX = 


( T( 66 )+T( 51 ) )/2 . 0 




TY = 


( T ( 67 ) -T ( 52 ) ) * . 7 5 + 


T( 52 ) 


TZ = 


( T( 67 ) +T( 52 ) )/2.0 




TAA = 


( T ( 6 8 ) -T ( 53) )*.75 + 


T( 53 ) 


TAB = 


( T ( 68 ) +T ( 5 3 ) )/2 . 0 




TAC » 


( T( 69 ) -T( 54 ) ) * . 75 + 


T ( 54 ) 


TAD = 


( T( 69 ) +T( 54 ) )/2 . 0 




TAE = 


( T( 69 ) +T( 70 ) )/2 . 0 




TAH = 


( T( 54 ) +T( 55 ) )/2 . 0 




TAF = 


(TAE-TAH) * .75 + TAH 




TAG = 


( TAE+TAH ) /2 . 0 




TAI « 


( T( 69 )-T( 70 ) ) * . 7 5 + 


T ( 70) 


TAL » 


( T ( 5 4 ) -T ( 55) ) * . 7 5 + 


T ( 54 ) 


TAJ « 


(TAI -TAL) *. 75 + TAL 




TAK = 


( TAI+TAL ) /2 . 0 




TAM = 


( T ( 70 ) -T ( 55 ) ) * . 7 5 + 


T( 55 ) 


TAN = 


(T(70)+T( 55) )/2 . 0 




TAO « 


( T ( 7 0 ) — T (71) ) * . 7 5 + 


T( 71 ) 
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TAR = ( T( 55 ) -T( 56 ) ) * . 75 + T(55) 

TAP - ( TAO-TAR ) * . 7 5 + TAR 

TAQ - ( TAO+TAR ) /2 . 0 

TAS - ( T( 71 )-T( 56 ) ) * . 75 + T ( 56 ) 

TAT = ( T ( 71 ) +T ( 56) )/2.0 

TAU - ( T( 75 ) -T( 60 ) )*.75 + T ( 60 ) 

TAV - ( T( 75 ) +T( 60 ) )/2 . 0 

TAW = ( T( 35 )-T( 34 ) ) * . 75 + T(34) 

TAX = (T( 36)-T( 35) )*.75 + T(35) 

TAY - ( T( 36 ) +T( 35 ) )/2.0 

TAZ - ( T ( 39 ) +T( 40) )/2 . 0 

TBA = ( T( 39 ) -T( 40 ) ) * . 75 + T(40) 
TBB = ( T( 40 )-T( 41 ) ) * .75 + T(41) 
TBC - (T( 5)-T( 4 ) ) * .75 + T( 4) 

TBD = ( T ( 6 ) — T ( 5) ) * . 7 5 + T ( 5 ) 

TBE = ( T ( 6 ) +T ( 5) )/2 . 0 

TBF = ( T( 9 )+T( 10 ) )/2 . 0 

TBG = (T(9)-T(10))*.75 + T ( 1 0 ) 

TBH => ( T( 10)-T( 11 ) ) * . 75 + T ( 1 1 ) 



:alculates average temperatures for rectangular element FEM program 
veld. for (chip and solder) 

TEMPSl = (T(95)+T(85)+T(96)+T(86) )/4.0 
DO 100 I - 1,115 

TEMPT ( I ) = TEMPSl 
TEMPB(I) = TEMPSl 
CONTINUE 

:alculates average temperatures for rectangular element FEM program 
veld. for (substrate) 



TEMPT( 116 ) 


= 


( T( 61 ) +T( 61 ) )/2 . 0 


TEMPB( 116 ) 


= 


( TA+TC ) /2 . 0 


TEMPT( 117 ) 


= 


( T( 62 )+T( 63 ) )/2 . 0 


TEMPB( 117 ) 


= 


( TC+TD ) /2 . 0 


TEMPT( 118 ) 


= 


( T( 63 ) +T( 64 ) )/2 . 0 


TEMPB ( 118 ) 




( TE+TG ) /2 . 0 


TEMPT( 119 ) 


= 


( T ( 6 4 ) +TI )/2 . 0 


TEMPB( 119 ) 


= 


( TG+TJ )/2 . 0 


TEMPT ( 120 ) 


= 


( T( 65 ) +TI )/2 . 0 


TEMPB( 120 ) 


= 


( TJ+TM ) /2 . 0 


TEMPT( 121 ) 


=s 


( T ( 6 5 ) +TO )/2 . 0 


TEMPB( 121 ) 


= 


( TM+TP )/2 . 0 


TEMPT( 122 ) 


= 


( TO+TS ) /2 . 0 


TEMPB( 122 ) 


= 


( TP+TT ) /2 . 0 


TEMPT ( 123 ) 


= 


( T ( 66 ) +TS )/2 . 0 


TEMPB( 123 ) 


= 


( TT+TW ) /2 . 0 


TEMPT( 124 ) 


= 


( T ( 66 ) +T ( 67 ) )/2 . 0 


TEMPB( 124 ) 


= 


( TW+TY ) /2 . 0 



91 



TEMPT ( 125 ) 




( T( 67 ) +T( 68 ) )/2.0 


TEMPB ( 125 ) 


=* 


( TY+TAA ) /2 . 0 


TEMPT ( 126 ) 


3 


( T( 68 ) +T( 69 ) )/2.0 


TEMPB ( 126 ) 


= 


(taa+tac)/2.0 


TEMPT ( 127 ) 


<1 


( T( 69 ) +TAE )/2 . 0 


TEMPB ( 127 ) 


- 


( TAC+TAF ) /2 . 0 


TEMPT( 128 ) 


3 


( TAE+TAI )/2 . 0 


TEMPB ( 128 ) 


= 


( TAF+TAJ )/2 . 0 


TEMPT( 129 ) 


= 


( T( 70 )+TAI )/2 . 0 


TEMPB( 129 ) 


= 


( TAJ + TAM ) /2 . 0 


TEMPT( 130 ) 


= 


( T ( 7 0 ) +TAO)/2 . 0 


TEMPB( 130 ) 


3 


( TAM+TAP ) /2 . 0 


TEMPT ( 131 ) 


= 


( T ( 71 ) +TAO )/2 . 0 


TEMPB( 131 ) 


= 


( TAP+TAS )/2 . 0 


TEMPT ( 132 ) 


- 


( T ( 7 1 ) +T (75) )/2.0 


TEMPB( 132 ) 


=s 


( TAS+TAU)/2 . 0 


TEMPT ( 133 ) 


3 


( TA+TC )/2 . 0 


TEMPB ( 133 ) 


= 


( TB+TD ) /2 . 0 


TEMPT( 134 ) 


3 


( TC+TE ) /2 . 0 


TEMPB( 134 ) 


= 


( TD+TF )/2 . 0 


TEMPT ( 135 ) 




( TE+TG ) /2 . 0 


TEMPB ( 135 ) 


= 


( TF+TH )/2 . 0 


TEMPT( 136 ) 


= 


( TG+TJ ) /2 . 0 


TEMPB ( 136 ) 


= 


( TH+TK ) /2 . 0 


TEMPT ( 137 ) 


3 


( TJ + TM ) /2 . 0 


TEMPB ( 137 ) 


= 


( TK+TN ) /2 . 0 


TEMPT ( 138 ) 


3 


( TM+TP ) /2 . 0 


TEMPB ( 138 ) 


3 


( TN+TQ ) /2 . 0 


TEMPT ( 139 ) 


3 


( TP+TT ) /2 . 0 


TEMPB ( 139 ) 


rs 


( TQ+TU ) /2 . 0 


TEMPT( 140 ) 


3 


( TT+TW ) /2 . 0 


TEMPB( 140 ) 


3 


( TU+TX ) /2 . 0 


TEMPT( 141 ) 


3 


( TW+TY ) /2 . 0 


TEMPB ( 141 ) 


3 


( TX+TZ )/2 . 0 


TEMPT ( 142 ) 


3 


( TY+TAA )/2 . 0 


TEMPB ( 142 ) 


= 


( TZ+TAB )/2 . 0 


TEMPT ( 143 ) 


3 


(TAA+TAC)/2.0 


TEMPB ( 143 ) 


= 


( TAB+TAD ) /2 . 0 


TEMPT( 144 ) 


3 


( TAC+TAF )/2 . 0 


TEMPB( 144 ) 


3 


( TAD+TAG ) /2 . 0 



TEMPT ( 145 ) 


= 


( TAF+TAJ )/2 . 0 


TEMPB ( 145 ) 


= 


( TAG+TAK ) /2 . 0 


TEMPT( 146 ) 


=r 


( TAJ+TAM )/2 . 0 


TEMPB ( 146 ) 


=* 


( TAK+TAN ) /2 . 0 


TEMPT( 147 ) 


a 


( TAM+TAP )/2 . 0 


TEMPB( 147 ) 


a 


(TAN+TAQ)/2.0 


TEMPT( 148) 


a 


( TAP+TAS ) /2 . 0 


TEMPB( 148 ) 


= 


( TAQ+TAT )/2 . 0 


TEMPT ( 149 ) 


= 


( TAS+TAU ) /2 . 0 


TEMPB ( 149 ) 


a 


( TAT+TAV )/2 . 0 


TEMPT( 150 ) 


= 


( TB+TD ) /2 . 0 


TEMPB( 150 ) 




( T( 46 )+T( 47 ) )/2 . 


TEMPT ( 151 ) 


= 


( TD+TF )/2 . 0 


TEMPB ( 151 ) 


= 


( T ( 4 7 ) +T ( 48) )/2 . 


TEMPT ( 152 ) 


a 


( TF+TH)/2 . 0 


TEMPB ( 152 ) 


= 


( T ( 4 8 ) +T( 49 ) )/2. 


TEMPT( 153 ) 


= 


( TH+TK )/ 2 . 0 


TEMPB ( 153 ) 


= 


( TL+T( 49 ) )/2 . 0 


TEMPT( 154 ) 


= 


( TK+TL )/2 . 0 


TEMPB( 154 ) 


= 


( TN+T( 50 ) )/2 . 0 


TEMPT( 155 ) 


= 


( TN+TQ ) /2 . 0 


TEMPB( 155) 


= 


( TR+T ( 50 ) )/2 . 0 


TEMPT( 156 ) 


a 


( TQ+TU ) /2 . 0 


TEMPB( 156 ) 


= 


( TR+TV ) /2 . 0 


TEMPT( 157 ) 


= 


( TU+TX ) /2 . 0 


TEMPB( 157 ) 


a 


( TV+T ( 51 ) )/2 . 0 


TEMPT( 158 ) 


= 


( TX+TZ )/2 . 0 


TEMPB( 158 ) 


a 


(T( 51 )+T( 52 ) )/2 . 


TEMPT ( 159 ) 


a 


( TZ+TAB ) /2 . 0 


TEMPB( 159 ) 


= 


( T ( 52 ) +T ( 53) )/2. 


TEMPT( 160 ) 


a 


( TAB+TAD )/2 . 0 


TEMPB( 160 ) 


a 


( T ( 5 3 ) +T ( 54) )/2 . 


TEMPT( 161 ) 


a 


( TAD+TAG)/2 . 0 


TEMPB ( 161 ) 


=* 


( T( 54 ) +TAH )/2 . 0 


TEMPT ( 162 ) 


a 


( TAG+TAK )/2 . 0 


TEMPB ( 162 ) 


= 


(TAH+TAL)/2 . 0 


TEMPT ( 163 ) 


a 


( TAK+TAN )/2 . 0 


TEMPB ( 163 ) 


a 


(TAL+T( 55) )/2.0 


TEMPT ( 164 ) 


a 


( TAN+TAQ ) /2 . 0 


TEMPB ( 164 ) 


= 


( TAR+T( 55 ) )/2 . 0 



0 

0 

0 

0 

0 

0 



TEMPT ( 165 ) 


= 


( TAQ+TAT ) /2 . 0 


TEMPB ( 165 ) 




( TAR+T ( 56 ) )/2 . 0 


TEMPT ( 166 ) 


= 


( TAT+TAV ) /2 . 0 


TEMPB ( 166 ) 


= 


( T ( 56 ) +T( 60 ) )/2.0 


TEMPT ( 167 ) 


a 


( T ( 46 )+T( 47 ) )/2 . 0 


TEMPB ( 167 ) 


= 


( T ( 3 1 ) +T ( 3 2 ) )/2.0 


TEMPT ( 168 ) 


= 


( T ( 4 7 ) +T ( 48) )/2.0 


TEMPB ( 168 ) 


= 


( T( 3 2 ) +T (33) )/2.0 


TEMPT ( 169 ) 


= 


( T( 48 ) +T( 49 ) )/2.0 


TEMPB ( 169 ) 


= 


( T ( 33 )+T( 34 ) )/2 . 0 


TEMPT ( 170 ) 


= 


( T ( 4 9 ) +TL ) / 2 . 0 


TEMPB( 170 ) 


= 


( T ( 34 )+TAW)/2 . 0 


TEMPT ( 171 ) 


= 


( T ( 50 )+TL)/2 .0 


TEMPB ( 171 ) 


= 


( T( 35 )+TAW)/2 . 0 


TEMPT ( 172 ) 


= 


( T ( 50 ) +TR ) /2 . 0 


TEMPB ( 172 ) 


= 


( T( 35 )+TAX)/2 . 0 


TEMPT ( 173 ) 


— 


( TR+TVJ/2 . 0 


TEMPB ( 173 ) 


= 


( TAX+TAY )/2 . 0 


TEMPT ( 174 ) 


= 


( T ( 51 )+TV)/2.0 


TEMPB ( 174 ) 


= 


( T ( 3 6 ) +TAY ) /2 .0 


TEMPT ( 175 ) 


= 


( T( 51 )+T( 52 ) )/2 . 0 


TEMPB ( 175 ) 


= 


( T ( 3 6 ) +T (37) )/2 . 0 


TEMPT ( 176 ) 


= 


( T( 52 ) +T( 53 ) )/2.0 


TEMPB ( 176 ) 


= 


( T ( 3 7 ) +T ( 38 ) )/2.0 


TEMPT ( 177 ) 


= 


( T( 53 )+T( 54 ) )/2 . 0 


TEMPB( 177 ) 


= 


( T ( 38 ) +T( 39 ) )/2.0 


TEMPT ( 178 ) 


= 


( T ( 54 ) +TAH ) /2 . 0 


TEMPB ( 178 ) 


- 


( T ( 39 )+TAZ )/ 2 . 0 


TEMPT ( 179 ) 


= 


( TAH+TAL )/2 . 0 


TEMPB ( 179 ) 


= 


( TAZ+TBA ) /2 . 0 


TEMPT ( 180 ) 


= 


( T ( 55)+TAL)/2.0 


TEMPB ( 180 ) 


= 


( T( 40 )+TBA)/2 . 0 


TEMPT( 181 ) 


a 


( T( 55 )+TAR)/2 . 0 


TEMPB ( 181 ) 


=* 


( T( 40 )+TBB )/2 . 0 


TEMPT( 182 ) 


= 


( T( 56 )+TAR)/2 . 0 


TEMPB ( 182 ) 


= 


( T( 41 )+TBB )/2 . 0 


TEMPT ( 183 ) 


= 


( T ( 56 ) +T ( 6 0 ) )/2.0 


TEMPB ( 183 ) 


= 


( T ( 4 1 ) +T (45) )/2.0 


TEMPT ( 184 ) 


a 


( T ( 3 1 ) +T ( 32 ) )/2.0 


TEMPB ( 184 ) 


= 


(T(l)+T(2) )/2.0 



TEMPT ( 185 ) 


= 


( T ( 3 2 ) +T ( 33) )/2.0 


TEMPB ( 185 ) 


= 


(T(2)+T(3))/2.0 


TEMPT ( 186 ) 


= 


(T( 33 )+T( 34 ) )/2.0 


TEMPB ( 186 ) 


* 


( T( 3 )+T( 4 ) )/2 . 0 


TEMPT ( 187 ) 


s 


( T ( 3 4 )+TAW)/2 . 0 


TEMPB( 187 ) 


- 


( T ( 4 ) +TBC )/ 2 . 0 


TEMPT ( 188 ) 


= 


( T ( 3 5 ) +TAW ) /2 . 0 


TEMPB( 188 ) 


= 


( T( 5 ) +TBC )/2 . 0 


TEMPT ( 189 ) 


= 


(T( 35 )+TAX)/2 . 0 


TEMPB ( 189 ) 


= 


( T( 5 ) +TBD )/2 . 0 


TEMPT ( 190 ) 


= 


( TAX+TAY ) /2 . 0 


TEMPB ( 190 ) 


= 


( TBD+TBE )/2 . 0 


TEMPT( 191 ) 


— 


( T ( 36 ) +TAY )/2 . 0 


TEMPB( 191 ) 




( T ( 6 ) +TBE )/2 . 0 


TEMPT ( 192 ) 


= 


( T ( 3 6 )+T( 37 ) )/2.0 


TEMPB( 192 ) 


= 


(T(6)+T(7) )/2.0 


TEMPT ( 193 ) 


= 


( T( 37 ) +T( 38 ) )/2 . 0 


TEMPB( 193 ) 


* 


( T ( 7 ) +T ( 8 ) )/2 . 0 


TEMPT( 194 ) 


= 


( T( 38 ) +T( 39 ) )/2 . 0 


TEMPB ( 194 ) 


= 


( T ( 8 ) +T ( 9 ) )/2 . 0 


TEMPT ( 195 ) 


= 


( T ( 39 ) +TAZ )/2 . 0 


TEMPB( 195 ) 


=r 


( T ( 9 ) +TBF )/2 . 0 


TEMPT( 196 ) 


= 


( TAZ+TBA)/2 . 0 


TEMPB ( 196 ) 


= 


( TBF+TBG ) /2 ♦ 0 


TEMPT ( 197 ) 


= 


( T ( 40 )+TBA)/2 . 0 


TEMPB ( 197 ) 


= 


( T( 10 )+TBG)/2 . 0 


TEMPT( 198 ) 


s 


( T ( 4 0 ) +TBB ) /2 . 0 


TEMPB( 198 ) 


= 


( T ( 1 0 ) +TBH )/2 . 0 


TEMPT( 199 ) 


= 


( T ( 4 1 ) +TBB ) /2 . 0 


TEMPB ( 199 ) 




( T( 11 )+TBH )/2 . 0 


TEMPT ( 200 ) 


= 


( T ( 4 1 ) +T( 4 5 ) )/2.0 


TEMPB( 200 ) 


= 


(T( 11 )+T( 15 ) )/ 2 . 0 


‘use if printout 


is desired 


DO 900 I - 1 


,200 


WRITE ( 13 , 


*) 


TEMPT ( I ) , ' ' , TEMPB 


3 CONTINUE 







END 
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